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ABSTRACT 


This study considers statistical modelling of a special 
type of images, called textures. Both causal and noncausal 
neighbourhoods are considered. We consider the 2D time 
series models for textures and the parameter estimation 
schemes for such models* We also discuss the frequency 
domain characterization of textures and application of PN 
arrays for this simulations. Among the non-causal models, 
we present two important ones| the simultaneous auto- 
regressive (SAR) and conditional Markov (CM), Each of these 
models is characterized by a set of neighbours, a set of 
coefficients and a noise sequence of specified properties. 

To avoid excessive computation, we use an approximate maximum 
likelihood scheme for estimating the parameters. Asympto- 
tically consistent decision rules are given for choosing 
an appropiate SAR or CM model. Model studies have been 
substantiated using synthetic textures. An image segmen- 
tation scheme based on SAR model is proposed and studied 
in detail. The programs generated for this study are also 
useful for images other than textures. 



CHAPTER 1 


INTRODUCTION 

The subject of image modelling involves the constru- 
ction of models or procedures for the specification of 
images. These models serve a dual purpose that they can 
describe images that are observed and also can serve to 
generate synthetic images from the model parameters. We 
will be concerned with the models of a specific type of 
images, the class of models for representing textures in 
image. Proper modelling of image has become important 
because, the efficient design of computer algorithms fox 
image analysis, processing and synthesis can only be done 
using the framework of an image model. 

There are four important areas of image processing in 
which texture plays an important role: classification [1,2], 
image segmentation [3], realism in computer graphics [4] 
and image encoding [5], Besides being an intrinsic feature 
of realistic objects, texture also gives important information 
on the depth and orientation of an object [6], Julesz [7] 
considers the problem of generation of familiar textures, 
an important one from both practical and theoretical view- 
points, Understanding texture is also an essential part 
of understanding of human vision [8], These considerations 



2 


have led to an increased activity in the area of texture 
analysis and synthesis. 

There is no universally accepted definition for 
textures. Part of the difficulty in giving a definition of 
texture is the extremely large number of attributes of 
texture that we would like to include in the definition. We 
consider a texture to be stochastic, possibly periodic, 2D 
image field. The most important attributes of texture 
are coarseness, contrast,^ directionality, line linkness, 
regularity and roughness. 

Most texture researches concerning texture modelling 
are characterized by the underlying assumptions made about 
the texture formation process. There are two major assump- 
tions, and the choice of these assumptions depends primarily 
on the type of texture to be considered in the study. The 
first assumption which is called the placement rule viewpoint, 
considers a texture to be composed of primitives. These 
primitives may be of random or of deterministic shape, such 
circles, polygons or even dot patterns. Macrotextures 
have large primitives, whereas micro textures are composed of 
small primitives. The texture image is formed from the 
primitives by placement rules which specify hov/ the primi- 
tives are oriented, both on the image field and with respect 
to each other. Examples of such textures include tiling 



3 


of the plane, cellular structures and a picture of brick 
wall. The second viewpoint regarding texture generation 
processes involves the stochastic assumption. The places 
ment rule principle for textures may include a random 
aspect in the stochastic point of view, however, we may 
consider that the texture is a sample from a probability 
distribution on the image space. The image space is 
usually a square lattice and the value at each grid point 
in the lattice is a random variable (non-negative). Textures 
such aS sand, grass and water are hot appropiately . described 
by a placement rule model. The key feature of these images 
is that the primitives are random in shape and cannot be 
easily described, 

1.1 TEXTURE MODELS 

By the model of a texture, we mean a mathematical 
process which creates or describes the textured image. The 
primary object of texture nradelling is the description 
of real textures* A secondary goal of texture modelling is 
classification of textures. The numerical parameters of the 
model can be used as features to classify the textures. 
Texture modelling is an important concept in the sense that 
it gives the inherent features of the image and hence helps 
further processing. A survey of the well known and commonly 
used textures model schemes has been given by Ahuja [9], The 
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problem of characterization of stochastic texture is basically 
a system identification problem where the system is a two 
dimensional one, 

1,2 THE PROBLe^ AND THE APPROACH 

In this thesis, we consider different approaches for 
characterizing a stochastic texture. Among the various 
methods of modelling stochastic texture, we only consider 
the random field (RF) models as they are mathematically 
tractable. In RF model of texture, we assume that the texture 
is generated by a linear system [lO] when the system is 
excited by white noise. The different pixel values of the 
texture follow a probability distribution depending on the 
system and the probability distribution of the input. This 
problem of system identification can be tackled in two 
ways 

i) Spatial Domain Approach 
ii) Transform Domain Approach 

In spatial domain approach, we represent the input-output 
relationship by a difference equation. Based on the depen- 
dency of any pixel on other pixels in a given region, we 
define a spatial model to be either causal or noncausal 
and they are shown in Fig, 1,1, 
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Fig. 1.1; Spatial Neighbourhood for images 


In case of noncausal model, we assume that the value 
of a pixel depends on all the neighbouring pixel values in a 
given region. In spatial domain characterization, the 
problem is to estimate the coefficient or the parameters 
of the governing difference equation. Different statistical 
schemes are generally used to estimate the parameters. For 
noncausal model, we assume the image to be torus, finite 
lattice image and depending on the input sequence, those 
noncausal models are of two type 

i) Simultaneous Auto Regressive (SAR) 
ii) Conditional Markov iOA) 

In SAR model, we assume that the input to the system is 
uncorrelated but it is correlated one in case of CM model. 
These types of noncausal models are very useful in texture 
analysis due to their realistic assumption about the 
neighbourhood. 
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In the transform domain characterization > texture is 
vievi/ed as the output of a homogeneous spatial filter excited 
by white, not necessarily Gaussian noise. Different stati- 
stical measures give the amplitude and phase response of the 
filter. In identifying the system, this method is not as 
good as the parameter estimation scheme but this is a 
computationally efficient scheme due to fast transform 
algorithms, 

1,3 LITERATURE SURVEY 

The random field model of textures has been studied 
by many researchers in the recent years. Depending on the 
type of image, the random field model has two subclasses 

i) Global models 
ii) Local models 

Global models treat an entire image as the realization of 
a random field. Hunt [11,12] has proposed a nonstationary 
gaussian RF model for textures* He has shown the appro- 
piateness of this model by subtracting the local ensemble 
average from each image point and showing that the resulting 
image fits a stationary Gaussian model. In a recent paper, 
Hunt [13], discusses three kinds of nonstationarities: 

Case 1: Monstationary mean, nonstationary autocorrelation. 
Case 2: Monstationary mean, stationary autocorrelation. 
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Case 3; Stationary mean, nonstationary autocorrelation# 

The nonstationarity in the image mean is described 
by the array of individual means taken over a specified 
neighbourhood about each image point. For the autocorrelation 
function, the breakdown of stationarity is related to the 
way the correlation function changes over the image. Hunt 
uses three attributes of the correlation function 

to describe its spatial dependence: 
its energy, its width, and its shape, Trussel and Kruger 
[14] claim that the Laplacian density function is better 
suited for modelling high pass filtered imagery than the 
Gaussian function. They have also shown that the basic 
assumption which allows the Gaussian model to be used for 
image restoration purposes are still valid under a 
Laplacian model, 

Nahi and Jananshahi [15] suggest modelling the texture 
by background and foreground statistical processes. They 
assume that the two portions are statistically independent 
random processes with known (or estimated) first two 
moments , 

Recursive solutions based on differential equations 
are common in one-dimensional signal processing and have 
been generalized to two dimensions, Jain [l6] investigates 
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the applicability of three kinds of random fields to the 
image modelling problem, each characterized by a different 
class of partial differential equations (PDEs): hyperbolic, 
parabolic and elliptic, A digital shape is defined by 
a finite difference approximation of a PDE, The class of 
hyperbolic PDE are shown to provide more general causal 
models than autoregressive moving average models. For a 
given spectral density function, parabolic PDEs can provide 
causal, semi-causal and noncausal representations. Elliptic 
PDEs provide noncausal models that represent two dimensional 
discrete Markov field and can be used to model both 
isotropic and anisotropic imagery, Jain argues that PDE 
models is based on a well-established mathematical theory 
and, further, that there exists a considerable body of 
computer software for numerical solutions. System iden- 
tification techniques may be used for choosing the PDE 
model for a given class of image, 

Pratt and Faugeras [17] and Gagalowicz [18] view 
texture as the output of a homogeneous spatial filter, 
excited by white, not necessarily Gaussian noise, A 
texture is characterised by its mean, the histogram of 
white noise and the transfer function of the filter. For 
a given texture, the model parameters are obtained as 
follows 
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i) The mean is readily estimated from the texture, 

ii) The autocorrelation function is computed to 

determine the amplitude response of the transfer 
function. 

iii) Higher order moments are computed to determine the 
phase response of the filter. 

Inverse filtering yields the white noise image and 
hence its histogram and probability distribution. The 
inverse filtering or decorrelation may be done by simple 
operators. For example, for a first order Markov field 
decorrelation may be done by Laplacian operator [19], The 
whitened field estimate of the independent identically 
distributed noise process will identify the spatial process 
in terms of autocorrelation function. 

If neither knov;ledge nor hypotheses about the type 
of random field are available, a class of local models i» 
used. These models assume relationships among gray levels 
of pixels in small neighbourhoods. The modelling process 
consists of choosing the relationship and evaluating its 
parameters. Two basic categories of local models arise 
from the joint and conditional probability formulations 
of variations in a neighbourhood proposed by Whittle [20] 
and Bartlett [2l], Besag [22], however described a 
number of difficulties with these approaches. He has shown 
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that the constraints on the conditional probability 
structure are so severe that they actually control the 
particular models* Nevertheless, the conditional proba- 
bility approach serves as the basis for a commonly used 
class of models, called Markov image models. Different 
choices for the set of a pixel’s neighbours give rise to 
different Markov models. The specification of the proba- 
bility distribution of a pixel’s gray level, given the 
gray level of its neighbours, defines a strict sense 
Markov field representation. Two of the more interesting 
models of this class have been discussed by Besag [23] and 
they are called autobinary and autonormal schemes, 

A close relative of the automodel is the simult- 
aneous autoregressive model( S/^Jl).KayshaP [24] notes that 
the symmetry requirements of the SAR are automatically 
fulfilled without the need to place prior restrictions 
on them, 

Kashyap [25] uses circulant (block) matrices to 
define SAR models for an infinite array (array obtained 
by repeating a given finite image). He has also given 
the scheme for estimating the parameters of the model 
and decision rules for the choice of optimum neighbours, 
Kashyap also discusses all these schemes for a doubly 
periodic array defined over a finite lattice. 
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The local models discussed so far relate the gray 
level of a pixel to the gray level of its neighbours, 
using a conditional probability formulation# Through this 
approach has got considerable attention, its primary 
difficulty, however, is the high dimensionality of the 
joint probability densities, even for small neighbourhoods, 
making parameter estimation complex and cumbersome# 

Rosenfeld and Troy [26] and Haralick [27] suggest 
the use of two dimensional spatial dependence of gray levels 
for fixed distances and angular perations, Haralick and 
numerous other investigators apply features derived from 
the co-occurence matrix to various texture classification 
and discrimination problems. 

The entire area of syntactic image model is not 
covered here as it is outside the scope of this thesis, 

1.4 OUTLINE OF THE THESIS 

After a brief introduction in Chapter 1, in Chapter 2 
we have studied the 2D time series model of textures. 
Synthesis and parameter estimation schemes for different 
auto-regressive and moving averages processes have been 
discussed, A procedure for checking the adequacy of the 
estimated parameter is also given# 

In chapter 3, we have studied the transform domain 
characterization of texture# Instead of computer generated 
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random number, simulation is carried out by using non-binary 
pseudo-random array, being ideally white in the cyclic sense. 
The amplitude and phase response have been estimated for a 
separable spatial filter. 

Chapter 4 discusses noncausal models for a finite 
lattice image. These images are assumed to be toroidal* 

Both simultaneous autoregressive (SAR) and conditional 
Markov (CM) images have been studied. For parameter esti- 
mation, we have followed Kashyap* s method for obtaining an 
approximate maximum likelihood estimate and a decision rule 
for optimum neighbourhood has also been studied. We have 
also discussed the spectral density characterization for 
such models* 

In Chapter 5, we have proposed an application of 
spatial interaction model to texture segmentation. Simula- 
tions have been carried out for different segmenting curves 
and for all such cases the proposed scheme worked satis- 
factorily, Some of the disadvantages and error performance 
of this method have also been discussed. Segmentation 
scheme is studied only for SAR models. 

Chapter 6 sums up the contribution of the thesis with 
recommendation for further work along the same line. 



CHAPTER 2 

CAUSAL SPATIAL MODEL FOR TEXTURES 

2.1 INTRODUCTION 

Time series modelling of any observed set of data 
is an old concept and is used extensively in statistics as 
well as engineering. These models are also being used 
successfully in social sciences and economics. Their success- 
ful application to ID data has led researchers to appl.ly 
such, models to multidimensional data} in particular tg 2D 
spatial processes. In the area of image processing such 
models are currently being used in image modelling and sub- 
sequently, using these models, to tackle problems of image 
restoration, texture classification, image data compression, 
spectral analysis etc. [ 9 ]. 

In recent years, however, considerable efforts have 
been made to model a 2D data array (e.g, images) by models 
which will take care of the dependency of a point on its 
spatial neighbourhood rather than on its 'past* as the 
concept of * past and future* are not defined naturally in 
such cases. 

In this chapter we will study two dimensional time 
series models and their parameters estimation, but this 
2D analysis will be done only after a brief introduction to 
its one-dimensional counterpart. 
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2.2 REVIEW OF ID TIME SERIES MODEL 

Let a time series be represented by • . 

and let ^ 0 )^^ represent an independent identically distri- 
buted (HD) sequence of zero mean and variance . Then 

the time series may be modelled by 


y . = p + (0^ + «l^i-l+«2^i-.2+ • - 


[I + S a. to 


k=o 


k^'i-k 


( 2 . 1 ) 


with = 1 
o 


or y^ may bo related to to^i^ by the eqn. 


^i " ^l^i-1 ^2^i-2 ••• + Wj, 


2 Piyi-i + 
i=l J ^ J ^ 


( 2 . 2 ) 


Here the eqn. (2.1) represents a moving average (MA) process 
of infinite order whereas eqn, (2,2) represents an auto- 
regressive process of infinite order. From eqn, (2.1) we can 
always get a moving average process of finite order (say q) ■ 
as 



q 

Z 

k-o 


“k“i-k 




(2.3) 
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Here jx represnts the mean of the output process. 

Again from eqn. (2.2) we can get an auto-regressive 
process of finite order (say p) as 


y, =* 


j=l ^ ^ 


•+W . 


(2,4) 


Eqn, (2.3) and eqn. (2.4) will be denoted as MA(q) and AR(p) 
processes respectively. By combining these two, we get an 
auto-regressive-moving average process where 


vi = t* + Piyi-i + 

j=i 


-I I"' ™ ^ ^ *? 1/' 

1 -• •• J 1“J K i-K 


(2.5) 


Eqn. (2.5) is often written as AF{MA(p,q) process. 

Assuming a zero mean process, eqn, (2.5) gives 


y«"^a*iy • it*.. .■t a*y • ^ ^ •'^■b 1^4 ^ 

1 1 1—1 p 1— P oil 1—1 


+ • • • + b X . 
f ... -r q 


(2.6) 


Relationship (2.6) can be represented by an input-output 
LTI model given by Fig. 2.1. 









2D LTI 


System 

' 


> 




Fig. 2.1; A time series model of Texture, 
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Taking Z~txansfoim on both sides of eqn, (2«6) 
+ ...+ a z'“*^)Y(z)-(bQ+b^z'*^+ b 


where the Z- transform is defined as 

oo 


Y(z) = s y-z 


-1 


So the transfer function of the LTI system is 


G(z) = If} 


^ -k 
b + S b, z ^ 

° k=l ^ 

P 

1+ E a. z 
k=l ^ 


for one sided Z-transform. 


Case 1: 


Let us assume bj^ = 0 for k = l,2f,.,q, 


yields 


G(z) = 


P -k 
1 + S a, z 


k=l 


By taking inverse Z-transform 

( 1 + E a. z~^ ) Y(z) ;= b.XCz) 
k=l ^ ° 


and 






we get, 

(2.7) 

( 2 . 8 ) 

Eqn. (2.8) 

(2.9) 

( 2 . 10 ) 


( 2 . 11 ) 
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which is an auto-regressive process of order p. So it is 
seen that in case of an auto-regressive process, the system 
is an all-pole system. 


Case 2: 


If we assume aj^ = 0 for k = 1,2,,, ^p, equation (2.8) 
reduces to 


G(z) = + E 


( 2 . 12 ) 


From eqn, (2,12) we get 


Y(z) = + E b. z"^ 

° k=l ^ 


(2.13) 


and 


= E b, X. , 
1 1, ^ k i-k 

k=o 


(2.14) 


This, in comparison to eqn, (2,3) is a moving average 
process of order q. 

Again it is seen from eqn, (2.13) that, in case of 
moving average process, the system is an all zero system. 
From eqn, (2,8) we see that the system in case of AHMA(p,q) 
process is a pole-zero system. In this respect ABMA comes 
out to be most general one. 


2,3 PRINCIPLE OF P.ARAMETERS ESTIMATION 

The auto-covariance for the case of AR(p) process 

is. 


( 2 . 15 ) 
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where E[,] represents statistical expectation. Then from 
eqn. (2.4) we get,by taking expectation on both sides 

Yk = ^iYk«i+^2”^k-2 ■*■•••+ Pp'Yk-p (2.16) 

The autocorrelation coefficient is defined as 

fk = <2.17) 

Hence from (2.16) and (2.17) we get, 

fk ^lfk~l ^2f’k-2 ^p'Pk-p (2.18) 

For k = 1,2,3, .. ,p, from eqn. (2.18), we get a set of linear 
equations 

fi “ + ••• ^pPp-1 

fa = ^ 2 * ••• ^p?p-2 

•i 

fp " ^ifp-i'^^afp-a'^ •••■*■ ^P 

This is the well known Yule-Vtelker [28 ] equations for 
AR processes. VJe can always find a matrix-vector relation- 
ship for eqn. (2,19). There is a fast algorithm developed 
by Levinson to solve such a set of linear equation for 
estimation the parameters [29], 

If we assume . be the ith coefficient in a kth 

K X 

AR process, then Yule Walker equation becomes 
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( 2 . 20 ) 


The left hand side matrix in eqn. (2.20) is of Toeplitz nature. 


so any fast method can also be used for solving such equa- 


tionSN. In the 
one, i.e.,0j^j^ 


set of coefficients <0^4$ f the last 

L ^^>i=l,,.,k 

is of great importance and it is called 


partial autocorrelation, and in an AR process 0j^j^ gives an 
estimate of the order of the process. From eqn. (2.20) we 
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In fact, partial autocorrelation function of any 
process reveals an interesting property of the model. It 
can be shown that PAC (partial autocorrelation) function 
^kk AEl(p) process follows 



for k £ p 

ii 

0 

for k > p 


Hence the PAG of a pth order AR. process has a cutoff after 
lag p. 

The moving average process, in comparison with AR 
process# is more complicated to analyse. Let us take the 
MA(q) process of eqn. (2.3) for the sake of analysis 


yi = yi - p 


( 2 . 21 ) 

Again, for such MA(q) process, the covariance is defined as 


Yk = E [YiYj__k] 

/ Yq = E [yiYi ] 

= +... + ^ 

=3 D (0^ 

q ' . • • • 

o o o 

where =3 (l+oCj^ -hcc^ + . ..+ cx^ ) 

2 

and = E[Wjj^a)^]= input variance 


( 2 . 22 ) 


(2.23) 
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Since we have assumed that the input sequence is zero mean 
and II D > it follows that 


= 0 

and E[co.tA).3 = 0 fox 1 4 2 
^ J 

From eqn. (2.21) it follows that autocovariance at lag k 
is given by 


Ifk = “l“k+l+ ••• + “q-k“q) 

for k = 1,2, •• .,q 

=0 for k > q (2.24) 

So combining eqn, (2,22) and eqn. (2.24), we get, auto- 
correlation 



lx = °‘k+”l°k-H+----" °‘a-k°‘a ) 

Vo Dq 


= 0 for k > q 


for k=l,2,..,q 

(2.25) 


From eqn. (2.25), for k = 1,2,3,..,, we get a set of non- 
linear equations as follows 

fl°q ••• ^q-l“q "7 

hs = + “q-2“q V (2.26) 

D = 2a„ 1 

Iq q q J 
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These equations can be solved for getting the 
estimates of the parameters. It is evident that^ for a 
lower order MA proioe-sei - thO • solution may be obtainable but 
when order of MA process becomes large, the solution becomes 
less tractable. 


An autoregressive moving average (ARMA) process is 
much more complicated than both AR and MA process in. the 
sense that it contains two sets of parameters 
In ARMA process, the parameter estimation is again done after 
computing the correlation coefficients. In case of ARMA(p,q) 
process, these will be q autocorrelations^^^, 
whoso values depend directly on the choice of q moving 
average parameters a as well as p autoregressive parameters 
Also the p values '•^he 
necessary starting values for the difference equation 
C^(B) = 0 ¥k > q+1, which entirely determines the auto~ 

correlation at higher lags. If q-p <0, the whole auto- 
correlation function l/'j = 0,1,2,..., will consist 
of a mixture of damped exponentials and/or damped sine waves, 
whose nature is governed by the polynomial C^(B) and the 
starting values. If however q-p >_ 0, there will be a 
q-p+1 initial values ^^q-p do not follow this 

pattern. These facts are useful in identifying mixed 


processes 
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Since parameter estimation in ARMA processes is not a 
simple procedure. Box and Jonkins have given a 

chart (P-520) which will give the initial estimate of the 
parameters for a ARMA(l,l) process, 

2.4 2-D STATISTICAL MODEL [30] 

In the previous section we have studied the analysis 
of different time series models, viz,, AR, MA and ARMA, 

Those methods can be applied to 2D discrete images only 
when the image array is converted to a 1^^ vector by some 
scheme (like lexicographic ordering). If z^^ is representing 
a Particular pixel value in the (i,j) position of the lattice 
then, it can be expressed again as a weighted sum of the 
data and white noise sequence (IlD) within a particular 
region or neighbourhood, A portion of such a lattice is 
shown in Fig. 2,2, and all, the grid points are equally 
spaced. 
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Fig. 2.2: 2D Image Lattice, 
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Hence in 2 dimension^^ the moving average process 
is given by 


2 • * * to • -] to • T •. a(0 • ^ « *1 to - « 4*^ • ♦ # * 

i»J 0,1 1,0 i-l,g 1,1 1-1, g 

(2.27) 

If we define the horizontal shift operator by H, then, 


H*^z - . — z. . 
ifl i,J-n 


(2.28) 


Similarly a vertical shift operator is defined by V where 



2 . 

i-ni ,3 


(2.29) 


Two more parameters need to be defined for the sake of 
analysis they are horizontal difference operator and 
vertical difference operator defined as follows. 


Horizontal difference Operators 

Vh=l-H 

■ Vh ^i,j = 0-H)z.^j= 

and vortical difference operator; 

Vv = i-v 


V . 


1,J 




4"*2. 

1,3 1 - 


■1,3 


(2.30) 


(2.31) 


These operators are commutative 



25 


z. . = V” z. . 


v:v: v: 


If we define 


(H.V) = E r Vn 

rnf=o n=o » 




(2,32) 


then eqn» (2*27) reduces to 

z,_j ^0.0 <2-33) 

Since '^(HjV) represents the system whose output is the 2^ 
time seriesj^z^ hence it may be called the system function 
of the model. 

Again in two-dimension, the auto-covariance at a 
lag (k,£) is defined as 


^k.t= 


(2*34) 


From eqn. (2#33) we get. 


Vlc.^ = 


OO OO C50 oo 


s Z Z S tp 

rio-.q=o m=o 


^^^“i-p,j-q ‘^i+k-m,j+4.- t ) 
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This expression is nonzero only wrfien pero-k and q=n-J2^ and 
for that value 


E[w. . f 0). . ] = cr 


oo oo 


. . ^k,l " P,q T P+k,q+^ 


( 2 . 35 ) 


The autocovariance coefficient ^ s't laq (k»£) can also 
be derived from the model system function ^(H,V). Let us 
define a function P*(H,V) 


A <y=> oo 

P(H,V) = X Z Y. ^ 

k=--« 1=^ 


H 


£. ^ 


(2.36) 


Then from eqn. (2.35) we get, 

P(H.V) = 


= “<0^ ' ?P.q 

= qj^^‘^(H,V)'$(H“^, v"'^) (2.37) 

Again for the sake of analytical simplicity we study 
the estimation technique for 2D AR and MA process. 
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Case 1; 


AR(m,n) process: 
Here (^H,V) 


m n 


\ III X* 1 

where Cp(H,V) = 2 E Q ^ 

k=o J^o 




and Cf = 1 

O jO 


Case 2: 


MA(p,q) process: 


Here z. . = ©(H,V) w. . 


where ©(H,V) =2 Z A 

k=o 1=0 


O3 s 0 3 ♦ 


ARMA (tn,mj p,q) process: 

Here q)(H,V) z. . = ©(H,V) w. . 

Xf J XyJ 


(2,38) 


(2.39) 


(2.40) 


jLike one dimensional ARMA process, the two dimen- 
sional ARMA process is complicated to analyse and the 
estimation of parameters in this case is not a straight 
forward procedure. In the remainder of this chapter 

we will study the AR and MA process in two dimension and 

will give an empirical rule to estimate the parameters. Test 
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for stationarity is also given for a given discrete 
data set. 

Let us consider an AB.( 1,1 ) and MA(l,l) nradel which 
are the simplest non-trivial 2-D processes. 


Case 1; 


AR(1,1) Process: 


S(H,V) = to, . 

i 1,3 

Here,. ^(H,V) = 1 + ^+<^0,1 ” ^1,1 ” ^ 


The- au-to-co variance is 


^k,l " ^1,0 ^0,1 ^l^l^k-l,,(,~l 

similarly 

;J*k,f ^1,0 k-1,^1‘^^0,1 k,!-!*^^!,! k~l, |rl 

Since, we are considering AR(l,l) process, hence by putting 
pairwise value of (k,l) as (0,1), (1,0) and (.1,1) we get 
the corresponding 2D Yule-Walker equation as 




29 


The two dimensional correlation coefficient at lag (k,x) 
defined as 

?k,£ ^0,0 

The given matrix vector equation may be solved for 
the initial estimate of the parameter set. It is evident 
that dimension of the matrix and vectors depend on the order 
of process. 

Case 2 : 

MA(l,l) process* 


Here = @ (H,V) 

0(H,V) = 1 + .0^qV + 0q^H + 

Since the autocovariance of z.j at lag (k,i^ is given by 


^k,l- ^i-k,j-l 

from the above equation we get 

Yoo = ^ ®10 + ®01'*' 

^10 "" ^®10 ® 11 ® 01 ^ 

^01 ^ ^®01 ® 11 ^ 10 ^ 

^11 %1 *^ 0 )^ 
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So autocorrelation functions are 




^10 ^lAi 
1+ + 0^2 




= ^n .t ^ 1^10 

ni o ^ ^ 

®oI 


fu 


9 


J-L 


2 « 2 2 
1+ 9j_q + ©01 + ©11 


From these equations we can solve for the parameters ©iq* 
9^, and 9-,. One thing to be noticed here is that the 
complexity of the parameter estimations increases with the 
order of the model, but in case of AR process if the order 
of the model is increased then the method for solving them 
remain unchanged in the sense that we have to solve more 
number of linear equations and same Yule Walker equation 
can be applied* So given a set of discrete data Whether it 
is one .dimensional or two dimensional, it is a general trend 
in all fields (where time series model works) to search 
for a suitable AR model, 

2.5 BAPIRICAL EQUATI.OMS FOR AODEX. BUILDING _ 

In application of these methods, to a given image 
data, the ensemble averages will be replaced by sample 
averages for computational conveniences* 
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For an image 
we define 
1) Mean 


ge ^ ^ij J where its dimensions is n^^ x n 2 » 


m ^2 


... 2 ^ 2^ Z * * 

m, n_ . , . , 1 T 

1 2 1=1 j=l 


2) Autocovariance 


'fk.t 


m-k 


where (k,£.) is the lag in {i,j) direction. 


3) Variance 


^ 2 
C = Y, 
Z ''' 


00 


4) Autocorrelation f) = 

J Yq^O 

5) Partial Autocorrelation function 


a) In horizontal direction 
► 1 


l>kk=? 0 ,.’ for k = 1 


k-1 f) 

I 0 ,k ' A f(k-i).al°'(''-^> 


1 ~ 






2 f 3 ^ • 
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b) In vertical direction 


tkk = fo,: 


k = 1 


fk,o “ f^(k~l),jf(k--j),0 

(k-l),j j^j,o 


k-1 

1 "• Z C|5 
3=1 


k — 2} • • I 


where 




ct, 


tk-2),3 ”'|^(k-l),(k-l) ^^(k-2), (k-j-1) 

'V'j = 1,2,, ..(k-1) 


2,6 MODEL BUILDIN3 

Given a set of 2D data we can not say by looking at 
it what type of model will be best for it. The information 
about the probable model (whether AH, MA or ARMA) can be 
obtained from the partial autocorrelation function. Box and 
Jenkins have shown [28] that the partial autocorrelation 
function gives an idea about the nature of the model as well 
as its probable order. They have given a chart of the 
nature of partial autocorrelation which is given hero for 
AH, MA and AHMA, 

We have also discussed in the previous section that 
a data set can be modelled as a time series model if and only 
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if it is a stationary one, A data set is said to be 
stationary if its covariance matrix is positive definite. 

The simple way to check stationarity of a series is by 
the behaviour of its autocorrelation function. In a 
stationary sequence the autocorrelation decays uniformly 
without any notch. In case of 2D data we have to test the 
smoothness of the autocorrelation in horizontal and vertical 
directions. With these things in mind we can give successive 
steps for model building for a given 2D data* 

Step 1 ; Calculate the mean value 

Step 2 ; Calculate the autocovariances 

Step 3 ; Calculate the. variance of the data sample 

Step 4 ; Calculate the autocorrelations 

Step 5 : Check stationarity of the given data, A given data 
may not be stationary in vertical or horizontal directions. 
This behaviour will be obvious from the autocorrelations. 

In most of the cases, fortunately ^it has been found 
that the da.fforoncGs in the adjacent . pixel values in 
the nonstationary process remove the trend and make the 
process stationary. So , in an attempt ( in general)to 
remove nonstationarity , we take differences in horizontal 
and vertical direction of the given 2D data till a stationary 
data series is apparent. 
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Step 6 ; Take the difference in the direction in which 
the data is nonstationary* The information about the 
direction of non-stationardiy can' b® obtained, from, the 
partial autdcorelation function. 

Step 7 ; Determine the partial ACF in order to check the 
approximate order of the model and check stationarity* 

Step 8 ; Identify the model fix>m Table 2.1. 

Step 9 : Estimate the parameters. After identifying a 
tentative model we can use the sample autocorrelations, to 
obtain a preliminary estimate of the parameters. These . 
values may be used as an initial estimate for obtaining maxi- 
mum likelihood estimate. 

For computing initial estimate of parameters in case . 
of AR pro cess, we have to solve a set of linear equations 
given by Yule Walker equation. But in case of MA parameter 
estimation, we have to solve a set of nonlinear equations 
and in that case determination of preliminary estimates is 
less tractable. 

After identifying a tentative set of model parameters 
for a given pictorial data, we will obtain the best-estimate 
of the parameters. This can be done by checking the rando- 
mness of the residual. The concept of residual in time 
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series model is as follows. Suppose, for any process 
(whether AR or MA),the actual value of parameter vector is 
©, Given a IID sequence (preferably Gaussian), we can 
always generate an image ^ by using the IID sequence and 

A 

vector 0. Let us say © is the estimate of the parameters. 

With the Same IID sequence and parameter © if we generate 
—an- image ^ then the difference between x Y is call! 
called the residual. If the model is best one, then tho 
residual Y^Cy^. = X “ Y.) will be white noise sequence. This 
is what will be achieved in step 10. 

Step 10 ; Compute the residuals and compute its autocorre- 
lation functions. The autocorrelation of the residual will 
suggest the direction in which the model should be modified. 

Step 11 ; Make modifications in the model based on the 
obtained information, 

. ■ * t 

All these steps are also given as a flow chart in Fig,2.3 
2.7 SIMULATION RESULTS; 

For simulation purposes, the textures were first 
generated using specific models and then the parameters 
were estimated from such synthetic textures. 

For this purpose, we have only considered lower order 
AR and Ma processes. Programs are also developed for 
generation of synthetic textures and subsequent parameter 




Fig, 2,2: Model Building Flow-chart 
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estimation and they are given in Appendix I, Provisions 
for handling higher order models is incorporated only, for 
AR processes because the nonlinear nature of the autocorre- 
lation equations in M/k process cannot be predicted from 
the knowledge of lower order MA models. 

In this study of 2D time series models, we have 
generated several synthetic images of size 64x64 and 32x32. 
Some of such synthetic textures are shown in Fig, 2.4, In 
these cases we have used only AR{1,1) and MA(l,l) processes, 
Gaussian textures were only considered and for synthesis we 
have used IMSL (Intn’l Mathematical Subroutine Library) 
routine GGNOR for generating Gaussian IXD sequence. But for 
2D time series models we need 2D, IXD sequence or an array. 
Since the elements of the sequence are independent so any 
arbitrary arrangement of them does not change their inde- 
pendence, Thus one way of generating a MxM UD array is to 
o 

generate a M xl HD vector and then arrange them (lexicogra- 
phically or by any other means) to form a MxM array. But 
considering the space requirements, v/e took another approach. 
In this method, each row (or column) of such IXD array, was 
generated by using different seed value of GGNOR routine. 
From here onwards whenever such computer generated random 
number array with prescribed density is generated we have 
adopted the same scheme. Though we have mentioned in the 
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last section that for best estimate of parameters we have 
to take the maximum likelihood estimate, but it has been 
found that for image size of 64 x 64 the initial estimates 
are quite close to the actual parameter values which is 
explainable from statistical - point of view. Instead of 
computing the maximum likelihood estimate which is a slow 
process in computers we have checked the whiteness of the 
residuals. Here we will give in a tabular form the 
actual and estimated parameter set for AR(l,l) and MA(l,l) proce 
process for a 64x64 image. We have also given the plot 
of the autocorrelation function of the residuals for finite 
number of points for AR(l,l) and MA(l,l) processes, in 
Figs, 2,5 and 2,6, 





Table 2*1: Summary of the Properties of ARy MA and AHMA Process [30] 
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Table 2,2 

IMAGE SIZE : 64 X 64 
MODEL : AR(l,l) 

Input ; GAUSSIAN (0,1) 


Actual Parameter Estimated Autoco- Estimated para- 

Values rrelations meter Values 
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0.821 4 
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0.6924 
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0.7964 


= 0.3601 

^01 

^11 

0.1832 

^Ql= - 0.1216 

= 0.1522 

= 0.2612 
(Zf^^ = -0.1013 
= 0.1567 
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CHAPTER 3 


FREQUENCY DOMAIN CHARACIBRIZATION OF STOCHASTIC TEXTURES 
3.1 INTRODUCTION: 

In the last chapter, we have studied the spatial 
domain characterization of textures by standard two dimen- 
sional time series models. It was also stated that models 
based on either autoregressive or moving average processes 
were computationally efficient. In general, parameter 
estimation for ABMA process is a complicated procedure. 

In this chapter, we shall study another procedure for chara- 
cterizing the texture. This procedure represents the system 
in terms of amplitude and phase response. This approach 
may be thought of as a frequency domain equivalent of time 
series model. 


White 
no ise 


Homogeneous 

Spatial 

System 


Texture 


Fig. 3.1; A model for stochastic texture 

The Figure 3.1 shows the model of texture generation 
Like the two dimensional time series model, here also we 
assume that stochastic texture is generated by passing a 
white noise sequence through a linear homogeneous spatial 
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system. In our study, the 2D system is a 2D digital filter# 

In literature, the white noise sequence is generally taken 
as the computer generated random number sequence or IID 
sequence, having some specified probability distribution. 

But, in tl^e frequency domain characterization of textures 
we use a 2D pseudo random array as the input noise array. 
Pseudo random sequences (generated from shift register 
sequences )have the property that they show nearly ideal 
autocorrelation function when defined in a cyclic fashion 
and they are also periodic. So, when such a doubly periodic 
array is passed through a linear system, the output also 
becomes periodic. So computation of amplitude and phase 
response can be done by FFT techniques. The useful concept 
for such system characterization is the moments of different 
orders. This method of texture modelling was tested here by 
generating a synthetic texture by such a process. For 
this simulation, nonbinary PN array was used as noise input. 

3,2 PN SEQUENCES AND ARRAYS; 

In this section, we will study the relevent properties 
of pseudo random sequences and its generation. It is found 
[31] that sequence, when generated as a maximal length 
sequence, shows some good properties. Such binary sequences, 
if they are generated by a m length shift register, will 
be periodic sequences of period 2'^--l, If we compute the 
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cyclic autocorrelation of such an n length sequence (ns=2®~l), 
the the autocorrelation follows the given relation 

•^(0) = 1 (3.1) 

j>(i) = ^ i , for 1 < i < n-l 

where '^(.) is the normalized autocorrelation function. Such 
a binary PN sequence approximates a white noise sequence, 
because, for a large length sequence, ^(i),-V'i 4 0 becomes 
negligibly small, which leads to a periodic delta function 
structure of autocorrelation. Such sequences have been used 
extensively in digital data scrambling, synchronisation etc., 
and a considerable body of literature exists [32,33,34 3. 

Recently,, several applications [ ] have called for two 

dimensional arrays whose 2D autocorrelation function is such 
that'^(0,0) = 1 and^(i,j) is small for (i,j) 4 (0,0). It is 
very easy to generate such arrays from a given PN sequences 
[31 3 • For generating such an array of dimension nj^xn 2 with 

Pco.o) = 1 

( 3 . 2 ) 

0 1 j < ^2-1 

(i,j) ^ (0,0) 

we need a sequence of length n, where n - nj^n 2 = 2”^-l, 
provided nj,and n 2 are relatively prime. For generating a 



45 


31x33 array (as in our case), we need a PN sequence of 
length 1023, Since the generated sequence should be maximal 
length sequence in a specific Galois field [ ]# we have 

to take the required primitive polynomial in that field. 
Primitive polynomial in a shift register sequence doteiy* 
mine's the governing difference equation. For the 

generation of a n-length shift register sequence (n=2™-l), 
we need m initial values to start with. So, given the 
initial values and a primitive polynomial in the specified 
field, we can always generate a maximal length sequence. 

The most important property of such maximal length sequence 
(MLS) is its autocorrelation function. The autocorrelation 
function '^(i) a real (or complex valued) sequence s^, 

®1* ®2' “* ®n-l length n is defined in cyclic sense as 


f 


(i) 


1 

n 


n-1 

x: s .s . ~ , 

j=o 


(3.3) 


1 — Of If Hh 2l f 


where bar represents complex conjugate and j(^ means 
summation modulo-n. This autocorrelation function i) 
is a periodic sequencej^i.e, (n+i) = ^(i). The auto- 
correlation function is shown . in Fig, 3.2. The auto- 
correlation function of a pseudo-random sequence of length 
n = 2’'''^-! is given by 
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1 

Fig* 3,2: Cyclic ACF function for a PN sequence of length n 
f (0) = X 

f (i) = - S i 1 i 1 2"-2 

Instead of binary PN sequence , having entries from 
0 and 1 only, we can also have a nonbinary sequence, taken 
from an alphabet of q symbols, for any q which is a prime 
or power of a prime. 

Let q be any prime power and let GF(q) be the Galois 
field with q elements. We need a polynomial g(x) of degree 
m with coefficients from GF(q) which i) is not the product of 
two such polynomials of lower degree, and ii) has a zero as 
a primitive element of GF(q'''^), say ^ , V/e call such a 
polynomial primitive over GF(q), In case of GF(4) (=GF(2 )) 
we have the elements ^0,1, ^ where m is a root of 
the primitive polynomial of degree 2 in GF(2). If we take 

such polynomial as x +x+l then, since w is root of itp 

o 3 

CO + ( 0+1 = 0 cind CO = 



1 . 



47 


Suppose, the primitive polynomial is 

g(x) =: .•* + (3.3) 

with h^^GF(q), h^ 9 ^ 0. The shift register specified by the 
above polynomial is as follows. In the Figure 3.3, the 



Fig, 3,3: General scheme for PN sequence generation 

boxes contain elements of GF(q), say The 

feedback path then forms qs 

® i+m “^ 111 - 1 ® i+m-l“^m-2® i+m-2’* • * * "^l® i+l'^o® i - 

(3.4) 

Eqn, (3,4) gives the recurrence relationship describing the 
output sequence. This is an infinite . sequence of period 
q'^-1, A segment of such sequence of length q'^-1 is called 
the pseudo random sequence over GF(q). As an example, if / 
q = 4, m = 2 and g(x) = x +x+a), we get the PN sequence 
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0 1 1 0)^ 1 0 to t») 1 w 0 0) 6)^ (3*5) 

2 

of length 4 -1 = 15* This sequence can be modified to 
a sequence 

011310221203323 (3,6) 

O 

by replacing w and w by 2 and 3 respectively. 

To find the autocorrelation function of a nonbinary 
PN sequence, we have to get a complex valued sequence 
from the real one in order to use eqn, (3,3). The unique 
transformation (from real to complex) which gives the 
autocorrelation as in binary sequence is as follows* 

If q is a prime or prime power then the conplex 
sequence a is obtained from a by replacing each Y^GF(g) 
by exp [iT-l j. After getting such complex sequence 

if we use the eqn, (3,3) then we get, 

j>(o) = i,-p(i) = i < i < 

3.3 PSEUDO-RANDOM ARRAYS; GENERATION: 

Though there arO umpteen number of ways of making 
a matrix out of a vector, only one method out of them is 
useful in making pseudo random arrays so that the array has 
properties similar to that of the generating sequence. 

In the context of pseudo random arrays it is needed that 
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this 2D array will have same autocorrelation structure as 
that of pseudo .random sequence. 

Firs-4 we will consider the binary sequence and 
corresponding array. A binary (modulo-2) PN array is 
constructed by folding PN sequence. Let us take number of 

k 1 krt k -1 

the form n=2 -1 such that n^^ = 2 -1 and n 2 = n/n^* 

Examples are 


n=15=2 -1 with K^=k2=2 nj^=3, n 2 = 5 

n=63=2^-l with Kj^=3, k2=2 = 7, 0^= 


Cyclic autocorrelation of such array (n^xn 2 ) is given by 

-^( 0 , 0 ) = 1 


The condition to be fulfilled is that n^^ and n 2 relati- 
vely prime. Here we will explain the construction procedure 


of such arrays. 


The starting point is a PN sequence a =(3^ taj.* • ‘^n-l^ 
obtained from a primitive polynomial of degree m=Kj^K 2 * We 
use a to fill up an nj^ x n 2 array by writing a down the 
main diagonal and continuing from the opposite side whenever 
an edge is reached. Figure 3.4 shows the construction of 
such array from a PN sequence when n=15^j^t=3 and n 2 =b. Thus 
the PN sequence 
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000100110101111 

produces the array 


0 

1 

1 

1 

1 

0 

0 

1 

1 

1 

0 

0 

1 

0 

0 

: 1 


(3.7) 


(3.8) 


Start at 
north-west 
corner and 
move south' 
east 



Continue from 
opposite side 
when an edge is 
reached 





^3 



®i 




^5 


^2 

1 



The complete 
array 


% 


^12 

CO 


^10 

^1 


®13 

^4 

®5 1 

^11 
j 

®2 

^8 

®14 


Fig. 3.4: Construction of a PN array 







These arrays are doubly periodic with period 
(horizontally) and rx^ (vertically). For such periodic 
array we can also define 2D cyclic autocorrelation function 
as 


R 


(k,j2) = 


^l”^ n2-l 


^ 1*^2 


Z Z x(i,j) X (i(+\ ^) 

i=o j=o W f^l W ^2 

(3.9) 

k = 0,l,..,,(nj^-l) 
and — 09 ly.,»,^( ) 


where x(i,j) represents the PN array and ^ is addition 
modulo n. 


For such a PN array, it is found that 
R(0,0) = 1 

and, R(k, ■i) = - ^ (0,0) 


(3,10) 


So, for large values of n^ and n 2 » the autocorrelation 
will be an approximate two dimensional delta function which 
will give an approximate flat power spectral density. 

Same procedure holds for generating non-binary sequences. 


3.4 SPATI/X FILTER MODEL 

In this section, we will try to model a texture as a 
homogeneous 

realization of a ‘ spatial process. Though it is 


52 


a standard procedure to create stochastic texture from 
computer generated random numbers, in this thesis we have 
studied textures generated from PN arrays. We make the 
assumption that the process is homogeneous and ergodic. 

Though homogeneity is not a valid assumption for real images, 
it may be a good one for texture type images. The spatial 
process or the texture is the output of a homogenous spatial 
filter,- Most of the time series models assume that the 
white noise sequence is a Gaussian one, but here we do not 
make any such assumption. In fact, roost of the real life 
textures are non-Gaussian, Since modelling of a random 
image means finding the system which will generate the 
image v\hen excited by white noise, we can represent the 
system in two v/ays, i) in sample domain and ii) in frequency 
domain. This particular model represents a texture (or 
homogenous image) in terms of the frequency domain character- 
ization of the system, 

3.5 TEXTURE ANALYSIS: 

In this section we will show that if a texture can be 
represented by the proposed model in Fig, 3,5, then it is 

^ Textufrb array 


Fig, 3,5: Proposed Model 


PN array 
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possible to compute all the filter parameters and then by 
inverse filtering we can get back the input noise 
parameters, 

a) Set of Texture Parameters: 

When a texture is generated or is supposed to be 
generated by the above model, then the resulting parameter 
sets for the textures are 

i) the texture mean 

ii) the histogram of the input white noise 

iii) the transfer function H(z) of the filter (or its 
impulse response). 

The second subset is equivalent to the set of different 
moments of the input IID sequence. These moments may be 
defined as 

Bn = nth order moment 
= E[b"(k)] 

where b(k) is the input noise. 

This set of moments and histogram of the noise are equivalent 
in that sense any one of them will completely characterize 
the noise sequence. 

3.6 parameter extraction [18] 

In this scheme of image modelling, the process of 
parameter extraction can be given in block form as in Fig,3,6. 
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Fig. 3.6: Block diagram for parameter estimation 

3,6,1 Amplitude Response Estimation: 

Let us assume t(kj^,k 2 ) be the pixel intensity at 
the point (kj^,k 2 )« Vt/e also assume that input array 
h{kj^fk 2 ) is doubly periodic with period Ml and M2 in 
horizontal and vertical directions respectively. Let us 
assume that h(m,n) represents the impulse (2D) response of 
the filter. So given h(m,n) and b(kjjk 2 )»WG can write the 
output t(kj^,k 2 )as the 2D convolution of h(m, n) and 
b(ki,k2)as 

t(ki,k 2 )= h(m,n) * * b(kj^,k 2 ) 
where * * represents the 2D convolution 
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ox. 


Nl- 

Z 


1 N2-1 

Z bCmj^,in2) hCkj^-m^^jk^-ni^) 


mi=o xa^o 


(3.11) 


Here we can also define 2-dimensional Z-transform of an 
array x(i,j) as 


oo oo . ^ 

= s Z x(i,j)z 2 (3,12) 

X=3 --.OO j=; — oo ^ 


So if we define TCzj^jZ^), H(z^, 22) and B( 2^,22) as the 2D 
Z- transforms of the array t(kj^,k2), ^^^1*^2^ b(kj^,k2) 

then from eqn. (3,11) we get, 

1(21,22 ) = H(2j^,22).B(2jl»Z2^ (3.13) 

Though the first order moment or mean of texture t(kj^,k2) 
is defined statistically as 

t(ki,k2) = E[t(kj^,k2)3 

it is estimated as 

— 1 ^ 2 ~^ 

t^(kl»k2^ = (3.14) 

The second order moment of the texture t(kj_,k2) gives 
the autocorrelation function array texture. 

The second order moment of t(kj^,k2) is defined again in 
statistical manner as 
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M, 




Using eqn, (3,11) we get 


N^-1 N^-l 


M, 


,(?.1»?2^“ S S b^(mj_,in2)h(k,-ra^,k2-n\2^ 

m^=o ni2=o 


h ( k 3 


(3.16) 


N^-1 . N^-l 

82 £ £ h(m^,in2)h(mj_+/(j^,in2-b^2) 

mj^=o m2=o 


(3.17) 


Eqn, (3.17) is true only if input array is a white one, 


Here 82 represents the 2nd order moment of the input 
noise. For time being, we can assume 82 = 1* 

If we take the 2D Z-transform of M 2 ( kA'>> we get 


M2(Zi,Z2) = 




^2^^1' ^1 ^ ^2 ^ 

(3.18) 


M2(Zj^>Z2^” H(z^,Z 2 ).H(zj^’"^, z”^) (by using eqn. (3.17)) 

(3.19) 

M2(2i,22)= 1 h(z^,Z2)|^ (3.20) 


So we see, in 2D case, the 2D Z-transform of the 
second order moment of the output of a system gives the 
amplitude square response of the system. 
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If the process is Gaussian then the 2nd order moment 
and first order moment contain all the statistical infor- 
mation of the process. By inverse filtering with a zero-phase 
filter we can get back the white noise. In that case, it 
is not necessary to compute the higher order moments. But 
in case the process is not Gaussian, it is possible to 

estimate the phase of the filter from higher order 

moments, 

3,6.2 Phase Computation: 

Phase information is contained in the non-Gaussian 
part of the higher order moments. In the simplest non- 
Gaussian case we can assume that the process is not Gaussian 
at the third order i.e, B^ ^ 0, We can define third order 
moment as 

It shows that the general filter structure gives a four 
dimensional third order moment function. In order to avoid 
long computation we describe the process of phase estimation 
by considering a special filter - called separable filter. 

If a 2D filter is a separable one, then it can be thought of 
as a combination of two one-dimensional filters. In case 
of separable filter, the, 2D impulse response h(m,n) can 
be written as 

h(m,n) = hj^(m)h 2 (n) 


(3.22) 
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or in frequency domain, 

HCz^jZ^) = Hj^(2j_)H2(z2 ) (3.23) 

where h^(m) <^i:^ H^( Zj^) 
and h^Cn) <^>H 2 (z 2 ) 

VJe will show here hov; to use the concept of third order 
moment in. the estimation of phase of a unidimensional filter. 
For this purpose let us consider the Fig* 3.7. 

t(n) 


Fig. 3.7: One dimensional model 

Let us define the third order moment of t(n) to be 
where 

%(^1»&2> = ECt(0)t(ij^)t(^) (3.24) 

Mo(iifL)= Bq h(ra)h(nn-lJ^)h(nH-22) (3.25) 

^ ^ m=o 

be 

Here we assume b(n) and hence t(n),to ^H-periodic. Taking 
Z- transform of 

'*-1 ^2 
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Using eqn, (3,25), we get 

* 

M3(z^,Z2) = B3H(z^)H(z2)H (3.27) 

jw, 

Let us assume ® and z^ = e , then from 

eqn. (3,27) we get. 


gl 3 (Wj^,W 2 )= + C((t 02 )-C((wj,+W 2 ) (mod 2it) 


(3.28) 


where 


A 

= phase of M 3 (zj^,Z 2 ) 


(^(a)) = phase of H(z^) 

^ = O or 1 depending on sign of Bo. 

h 

So, after computing the phase of M 2 (zj|^,Z 2 ),we have 
to find the phase of the filter H(z) i.e. 0(w) by using 
eqn. (3.. 28). Here we will give an iterative procedure to 
compute the phase. 

Let us assume that we want phase (il5(a)) of an array of 
N values, depending on the size of the filter. 


Let us define = §2^!^ 


and 


211 


“2 ~ N 


(■t-i) 


Hence from eqn. (3,28) we get 
C<o(k,^)= -rc+0(k )+«((-£) - gJ(k+£-l) (mod 211 :) 


(3,29) 
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The eqn, (3,29) is the digitized version of eqn. (3,28), 
Forj^= 2 in eqn, (3,29), we get, 

03 (k, 2 ) =11 + 0(k)+0(2)-af(k+l) (mod 2it) (3,30) 

C^3(]>»2) = 11 + 0(l)+C((2)-eJ(2) (3*31) 


0(1) = 03(1,2) 


assuming 




Similarly, 


1 W 

0(2) = i E 03 (k, 2 ) 

I:C~1 


(3.31) 


(3.33) 


We can rewrite eqn, (3,30) (with ^= l) as 

0(k+l) = 0(k) + 0(2) - 03(k,2) (3.34) 

Vi/e will use eqn, (3,29) to eqn, (3,34) for the phase 
estimation of the 2D filter. 


3.7 USE OF SEPARABLE FILTER TO COMPUTE PHASE 

Vi/e have already stated that if the linear filter is 
not a separable one then third order moment yields a foul - 
dimensional function and from such function phase compu- 
tation is a very time consuming process. For this reason, 
we have taken a separable filter. In such cases, the 
convolution also becomes very convenient one. 
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A filter h(i,j) is said to be a separable one if 
vje can write 


h(i,j) = h^(i)h^(j) 
in Z-domain H(zj^,Z 2 ) = H^( Z 2 ) 


(3.35) 

(3.36) 


If we say X(z^,Z 2 ) is the Z-transform of the input and 
Y(Zj^,Z 2 ) is the Z- transform of the output, then entire 
process can be shown as in Fig. 3.8. 



Fig. 3.8: A Separable filter structure 

Though we can get Y by 2D convolution of X with H , 
still we can utilise the separable nature of. the filter. 
In the Fig, 3,8 the filters Hq and are both ID filter. 
Let us assume X is a n^^ x array. Then the first filter 
Hq will take each column of X in succession and will modify 
them to get another intermediate array Z, . Each column 
of Z is the column filtered- version of X. This new array 
J^is fed to the next ID filter will take each row 

of Z' and it will modify the row to produce a row of the 
array Y, So Y is the row-filtered version of the array Z.' 
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This separable property of the system can be exploited in 
computing the phase. If represents the phase of 

the ID filter and represents the ID phase of Hq 

then phase of the total filter H is given by 

(3.37) 


So effectively the 2D phase of the system is computed by 
the help of its one-dimensional portions. 


3.8 SIMULATION RESULTS: 


In the simulation process we have taken a 31x33 doubly 
periodic pseudo random array and a 2D filter with transfer 


function given by 


2-2, 2-z^ 

Hiz^tz^) = (xi2^) ( 1-2 22^ 


(3.38) 


This filter is a constant magnitude (all pass) filter. It 

has the property of introducing no correlations. By 

“jw, -30)2 

replacing = e ^ and Z 2 = © we can compute the 
amplitude response of such filter on the unit circle. Then 


-jw. 


H(a)i,o)2) 


L2rS tl 


■J0)r 


.( 2 ~ e . IL .- 


(l-2e 


-JO), 


•) (l-2e 




(3,39) 



-Jto, 


Since, 

e = 

cos - 


“3^2 


and 

e = 

cos 0)n - 


(3.40) 



63 


simple trigonometric identities lead to 
lH(Wi,W2)i = 1. 


(3.41) 


This type of filter shows a pole-zero structure in its 
transfer function which in sample domain, is equivalent to an 
ARMA process. Since H(_z) = where N(z) and D(z) are 

polynomials in z such filter is an HR (infinite impulse 
response) filter. But impulse response of such filter is 
also a decaying one in both direction. The input v/hite 
noise array was taken to be 31x33 from GF(4), i,e, the 
elements of the array are from the set |^0,1,2,3'^ • Since 
total number of elements in this array is 1023 (=4 -l) we 
need a shift register of length 5, Vie also need a 
primitive polynomial of degree 5 in GF(4). The polynomial 
is taken to be g(x) 

g(x) = x'^ + X + w (3,42) 


where w is a root of the primitive polynomial of degree 2 in 

2 

GF(2). If we take that polynomial as x +x+l = 0 then 


p 

w + w + 1 = O 
and co^ = 1 


(3,43) 


So the required shift register is as shown in Fig. 3.9. 
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Fig. 3,9: The feedback connection for Eqn, (3,42). 

After generating the periodic sequence of length 1023 
we then used the scheme for converting a vector into a 
matrix to form a 31x33 array. This 2D array xCsay) was then 
filtered by the spatial system H(z2j22^ get the periodic 
output which Was the required texture. Though we have 
represented the filter in Z-domain by the transfer function 
H(Zj^,Z 2 )f the filtering scheme was developed in sampling 
domain by writing the corresponding difference equation. 

We also took advantage from the separable property of the 
filter. Here we refer to the Fig, 3,8, For the given 
system 


Hq { z ^ ^ 
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If we assume, 

«c‘Zl)= (3.44) 

then, by taking Z~ transform of eqn* (3,44), we get, 

y*(n)=0.5 x*(n)-x*(n-l)+0.5 Y*(n-l) (3.45) 

As mentioned earlier, Hq will take each colunaa of the 
input and produces a new column following eqn. (3.45). If 
z represents the 20 output of Hq then, 

z(if j)=0-£'x(i, j )-x( i, j-l)+0.5z( i, j-1) (3*46) 

Similar expression is obtained in the case of second filter 

y( i, j)=0.5z(i, j)-z(i-l, j) + 0.5y(i-l,j) (3.47) •• 

In using equation (3.46) and eqn. (3,47) it is to be 
remembered that the input is periodic and system is causal, 

i.e,, 

x(i,a) = 0, Yi < 

and 

y(i,J) =0, Vi < 0, j < o 

For computing the impulse response array, vje excited 
the system by a 2D delta function ^(i,j). For the compu- 
tation of 2D Z-transform (on the unit circle) of second 
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and third order moments of the output, we have used fast 
fourier transform routine FFTP (IMSL)* Fig* 3,3:0 through 
3,13 depict the 2D cyclic autocorrelation of the input, the 
actual and estimated amplitude and phase . response respe- 
ctively of the all pass system used in the study. In such 
cases we cannot take input noise array of arbitrary dimension* 
In order to maintain the idealness of autocorrelation, we 
have to arrange the elements of the sequence in a particular 
vjay to get the array* Considering all these facts, it seems 
that for such simulation purposes, computer generated random 
number are more convenient than shift register sequences. 

Such pseudo random sequences/arrays, can however still 
be used as a substitute for computer generated random 
number for many applications. The periodic structure and 
almost ideal autocorrelation of pseudo random numbers help 
us in this respect. 

Since input noise is a pseudo-random array, the 
generated texture should be named ‘pseudo random* texture. 








CHAPTER 4 


ANALYSIS AND SYNTHESIS OF TEXTURES BY f-JONC;jJS<:jL RF MODELS 


4,1 INTRODUCTION 

In the last two chapters, we have mentioned sample 
domain and frequency domain characterization of textures. 

In this chapter, we will study a new class of RF models 
where the concept of neighbourhood is a general one^ i,e,, 
it can also take of non-causal neighbourhood. These models 
may be called as non-causal spatial interaction model as 
described by Kashyap [lo ]. Such RF models represent an 
image as a distribution of gray levels over a finite or 
infinite lattice where each point is supposed to be equi- 
distant. 

Every such model is characterized by a set of 
neighbours, the corresponding coefficients and an independenti 
identically distributed random sequence with prespecified 
probability density. Typically an image is characterized 
by 2D scalar data, the gray level variations defined over 
a rectangular or square lattice. In such lattice structures, 
gray level at a particular lattice point has a statistical | 
dependency on the points inside a specified neighbourhood. 

There are two nonequivalent ways of specifying the 
interaction between pixels in a given lattice, leading to 
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two different desses of models, the simultaneous models 
[22 3 and conditional Markov models. Here we will 
discuss only the spatial interaction models over square 
lattice. These models dp not require that the underlying 
field is Gaussian. These models are therefore quite general, 
i.e,, a variety of image types can be generated by choosing 
different neighbour sets, IID sequences of different 
probability density functions and different parameter sets. 


4.2 NOTATIONS 


r = ^ ~ ^®1*®2^* Cn~n^dinates of grid points, 

r^ and s^ integers, y(s),= gray level or intensity at pixel 
position s, y(s) is a non-negetive, real number. 


SI 




= (Sj, , 32 ), 1 < Sj, ,S2 



=s set of grid points of a finite image of size MxM. 
= S = (Sj^jS^), - ~ < «» 

= col [y(l,l), y(l,2),..., y(M,M)3 


N = 
z = 



N 


vector of a finite MxM image 
a neighbour set, a subset of (0,0)^ 

(zj^,Z 2 )» a pair of operators, defined below 

If r,s^l2, zVs)^y(Sl+^l»S2+3r2) 

=: a pair of frequencies in both directions. 
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> i = frequency in ith direction 
expdPiT^) = (expYZi;\^^ 

A(z) = 1 - £ Q 

T'Sn ^ 

A^(r)=. A(z= e^27ir/M ^ 

¥(•) = D.F.T. of l^y(s), s ^ Jl_^ 

U(.) = D.F.T. of s ^J1 

R(*) = Auto covariance function 
S(A) - Spectral density 
p( . ) » Probability density 

The convention for labelling members of N are as follows; 

x(-l,0) North 

x(0,-l) V/est . x(0,l} East 

x{l,0) South 

For example if the set N is represented by the 
coordinate as given 

N = (-1,-1), (-1,0), (1,0), (1,-1), (1,0), (1,1) 

then pictorially, 
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4,3 SUMMARY OF DISCRETE SPATIAL INTERACTIOM MODELS 

In this chapter we will discuss only discrete spatial 
interaction models with gray levels either real values or 
integers. In all such models the input is a set of indepen- 
dent random variables u(s) with a common density|;)^. Such 
sequences are called independent identically distributed 
sequence (IID). The random field (RF) image y(.) is 
generated by transforming the input sequence in a suitable 
manner. In general, the system representation of such 
transformation may be of the following type « 



Fig. 4.1; A generalized model for RF Textures 

Firstly, we vvill discuss about T 2 end then with proper 
context we will discuss Tj, . The system T 2 is characterized 
by transfer function A(z) which is defined by a neighbour 
set N and postulated to be 

A(z) = 1 - £ 

r4N ^ 

The system Tj^ and T 2 may bo interpreted as a difference 
equation in sample domain or multiplication in transform 
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domain* The output of the first system Tj|^ may be corre- 
lated or uncorrelated depending on the system T^* Now if 
we Want to fit such a model for image data, then our main 
objective is to estimate the finite dimensional parameters 
the probability |)y( . ) of the input IID sequence. 
We can divide the generalized model (Fig, 4.1) in two ways: 


The infinite lattice random field (ILRF) models, 
where y(s) is defined over an infinitely extended grid or 
lattice, i.e., for such image y(i,j) - oo < (i,j) < «» . 


Ca se., 2 : The finite lattice random field (FLRF) model is that 
one whore y(s) is defined only over a finite region and in 
that casej’^= ^(i, j), l < i,j < representing MxM image* 

In general for RF models whether it is a infinite. ■ 
lattice or a finite lattice we can represent the image 
data y as 


y(s) 


a 


I ©u [y(s+ s, ) - 


Sj^€N 


s 


a] + VJu(s) 
(i,j) 

0 1 jh i f + 2 jr 


(4.2) 


Depending on the statistical nature of u(s) we may obtain 
various types of RF models such as SAR (Simultaneous Auto— 
Regressive), CM (Conditional Markov) etc. One advantage of 
such RF modes is that they incorporate geometrical and 
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statistical dependency of the pixel on its specified 
neighbourhood, not considered in many other models. By 
such models we can generate synthetic textures obeying 
a given model with known parameters by using a set of 
random numbers, 

4.4 INFINITE LATTICE JvDDEL 

Let us consider Fig, 4,1 as the general model. In 
IL models the transformation T 2 is interpreted as a difference 
equation whose coefficients are from A(z) and 

y(s) - £ 0 ^ y(s+r) = Vp\ 9 (s) (4,3) 

This equation may also be written as 

A(z) y(s) =l/^'\$(s), z =: (4.4) 

where 2 ^ are regarded as delay operators in the 2 directions 
and A{2) is given in equation (4,1), The different choices 
for Tj^ lead to different types of models. 

Simultaneous AR (SAR) Models : 

For a SM model defined over infinite lattice, 

we have 

Tj_ = 1 and v(s) = u(s) 

S imul ta noou s /JRMA • ( SARM<\ ) Mo dels ; 

For a S/iRMA model defined over infinite lattice, 
we have , 
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T 


i 


- B(z) 


EC*, are real 




I’ 


Equ ivalently. 


Y(s) - E y(s+x) 


*) 


x,N ‘ 

= V3[u(s) + E u(s+r) ] 
T r N ^ 


Separable Models ; 

Here A(z) and B(z) possess a factorization 

A(z) = 

B(z) = 


The neighbour sets characterizing A(z) and B(z) are located 
in a quadrant. 


Conditional Markov Model ; 

In conditional Markov (CM) models, the neighbour set 
H, used in characterizing the polynomial A and B is symmetric, 
i.o., rC,N=^-r^.N. Secondly the system is implicitly 
defined by the following relation. 

S^()^) = Spectral density of v(.) 

» A[expf*-1 2ic)^, exp f-l 271^2^ 


( 4 . 5 ) 
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Equxvalsntlyj the sequence v(*) has the following 
covariance function, 

Ry(s) = 1 if s = ( 0 , 0 ) 

~ — O if s^N (4,6} 

s 

= 0 otherwise 

For such models, the coefficients follows the 

relation 

0 ^ — 
s s 

V\/e can assume the relationship between u(.) and v(,) as 
follows 

= VaH) (4.7) 

To validate the assumption we can consider two mutually 
exclusive cases. 

Case i; There exists a polynomial Aj^Cz) so that 
A(z) = Aj_(z~^) 

and 

A,(z) = 1 - E , 0^ 2 ^ (4.8) 

^ rCNj 

In such a case T^^jin eqn. (4.7) is interpreted as a 

\ 

difference operator 

v(s) A^Cz) u(s) 


(4,9) 
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So in that case, the conditional model reduces to a simult- 
aneous model, 

assume that Aiz) has no factorization, so 
in that case eqn, (4,7) has no meaning. 

It is also important to note that a symmetric poly- 
nomial in a single variable is always factorable. 

The general nonfactoribility of a symmetric 2D polynomial 
AiZj^rZ^) is one of the key differences between the ID and 
2D theory, 

4.5 FINITE LATTICE MODELS 

These models were introduced by Kashyap [32 ] and 
they have many elegant properties. In FL models both the 
inputs and output (image) are assumed to be defined over 
a MxM grid. 

As we are interested in finite lattice image so in 
the subsequent sections v/e will discuss only FL model, 
however, an equivalence of these with IL models will also 
be mentioned. 

Neighbour Sets; Markov Definition and Stationarity 
Neighbour Sets ; A neighbour set N is any finite subset 
of the total set < « . For all such 

sets (0,0) Au A neighbour set is called symmetric if 


(i,j)€N ^(-i,-j)^N 
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The first order and second order neighbour sets 
are defined below. 

^ XXX 

X . X X . X 

^ XXX 

N: First order Ns Second order 

In the last two chapters we have analysed texture 
where neighbourhood was assumed to be causal* A causal 
neighbourhood has all its members in a quadrant 


N = (i,j), (0,0) N, i < 0, 0 < 0 


X 

X 


X 

X 


X 

X , 

X 


N: Causal 


N: Semi causal 


N is called as one-sided or unilateral if N is a finite 
subset of a half plane defined as 


Dofinition ; (Half Plane S'^): S"^ is recursively defined 


as follows 

i) SjiCs"*" , ^ 

iii) ( 0 , 0 ) ds* 
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From the above definitions it is seen that is not 
unique. Some examples of are 



The non-causal neighbourhood is a more general neighbourhood 
and their use is more realistic for discrete image data 
modelling. 

Markov Definitions ; 

ID Markov Sequences ; 

A unidimonsiondl sequence ^y(t), t = 1,2,.,.^ is 
said to be markov of order m if 

p(Y(t)ly( tjL),tj^<t)=:p(y(t)ly(t-l),...y(t-m)) (4.10) 

probability 

where p is the ^ density function. Here m is 

called the memory of the process. A ID process y(t) is 
called bilateral markov of order m if 

p(y(t)|all y(tj^), t^ t) 

= p(y(t)|y(t+l), y(t + 2),...,y(t+m)) 

(4,11) 
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In a similax inanner, we can define vector maxkov 
sequence* In such a case, every component of the vector 
is said to be projection of a maxkov process, 

2D Markov Definitions 1*351 

A 2D sequence y(s) is said to be bilateral maxkov 
with respect to a symmetric neighbourhood N if 

p(y(s)|all y(s»), s s') = p(y(s) jy(s+s< }, all s*^N 

(4.12) 

Comparing eqn. (4.11) and eqn. (4,12) we see that 
bilateral markov in 2D is just a generalization of bilateral 
maxkov in ID and in general bilateral marlq^v in 2D is 
known as 2D markov sequence. 

Stationarity and Isotropy ; 

A sequence y(s) is covariance stationary if 

Covariance [y(s), y(s+r)] = function of r only 
^ R(r ) = R(-r) 

where covariance [xy]A E[(X-E(x))(Y-E(y))] 

A sequence y(.) is weak stationary if it is 
covariance stationary and 

E[y(s)] ^ 

The sequence y(,) is called a strong stationary sequence if 
its probability density function obeys the following 


conditions. 
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p(y(s') = y(s") x^) = f(Xj^,x2,s'*-s'*) 

In this study by stationarity we will mean the 
weak stationarity in the covariance sense. A stationary 
sequence will be called isotropic if, 

R(s = (i,j)) = R(|il,)j|) = a(j,i) 

For an infinite lattice covariance stationary R.F, the 
spectral density S(^ is defined as 

S( =s I R(s} exp(V~l 2% 
sfcl^ 

oO 

- 2 K R(s,,s )exp[f-l27C(s > +s^ )] 

8^= -00 -CO ^ X i Z Z 

4,6 ANALYSIS OF FINITE LATTICE MODELS 

We have discussed in the previous sections some of 

the properties of infinite lattice RF(ILRF) models. The 

greatest impediment to the use of RF models defined over 

infinite lattices for analysing finite images is the 

extensive computation involved. Even the computation of 

2 2 

some of the basic quantities like M xM covariance matrix 
for an MxM image requires extensive computations. Conse- 
quently, most of the RF models, used in image processing 
use simplistic assumptions like isotropy, causal neighbour- 
hood etc,, so as to reduce the computational complexity. 
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Of course, there is considerable interest in the 
approximation methods. However, every approximation method 
seems to be tailored to a particular RF model and a 
particular application. 

The approximation of a covariance matrix by a 
circulant matrix or by block-circulant matrix is justified 
only by asymptotic considerations, i.e., the size of the 
image tends to infinity. These approximations are not 
suitable for synthesis of any image, obeying a given RF 
model. 

Here in this section, we will analyse the basic 
features of finite or lattice RF (FLRF) models for images. 

We will discuss the development of several familieis of 
RF models defined only over a finite square lattice, say 
Sit which is same as that of the given MxM image to be 
analysed. The basic idea behind the FLRF irodel is as 
follows: 

Using a given neighbourhood set N, characterizing the 
model or system, partition the given image (MxM) lattice 
into two mutually exclusive subsets XL| andJTIg. The 
boundary set Jig has every grid point which has at least 
one neighbour defined by N not in_X2>. .Ttj is the interior 
set, equals for generation 
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of y(s)> and y(s), s are different. Specifi- 
cally every pixel in follows the eqn. (4.2) l.e. 

they are same as that of ILRF model but each pixel in i\ « 
follows a modifiod equation as shown by eqn# ( 4 # 13 ) 

y(s) y(g^) +V^u(s) (4.13) 

® €AlQf S ^Jl , Sy^ ^ s 

Thus the equation for the image vector can be written as 

(4.14) 

where B(©) is M^xM^ and u is a M^xl vector of priraiti-^^ 
random numbers independent or dependent on the basis of 
SAR and CM models. By choosing appropiately, 

the matrix B(0) can have eigenvectors, either fourier 
vectors or discrete sine or cosine vectors as the case 
may be. The distinction between andjflg can be given 
as follows 

j~i.g » |^s~(i,j), s^/Uand (s+Sj^)-^ J\ for at least 
one number Sj^£: N ^ 

ill =* |s=:(i,j),s AJl and ( s+s^^) 

The eqn. (4.13) was chosen for the pixels y(s), 


such that 
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i) the equations axe closed, i,e, they do not involve 
any s’ where s» ^JTL 

ii) the overall set of equations are easy to analyse, 
SAR Models : 

On the basis of the above mentioned partition of the 
finite lattice into XI j g we can write the eqn, 

^4*2) as 

for y(s) obeys (4,2) 

but for s£Jl,g 

y(s) - a = « (s + k) - a]+Vp u(s) (4,1b) 

where 

y^(s+k) with s-(i,j) and k=(p,q) 

- y(s+(P,q)) if (s+(|>,q)) C. 

= yCCp+l~l)j^+li(q+j-l)^+l] (4.16) 

if s+(p,q) ^ Jl, 

This modified eqn, (4,16) is based on the assumption that 
the lattice is a torus one, i,e,, for such lattices the 
opposite edges are similar. This is a very useful conception 
used in this and subsequent chapters, 

Eqns, (4,16) and (4.2) give eqns. relating the 
image variable y(,) and the IID variable u(,). By using 
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the same relation as given in ec(n* (4.14) we see that for 

toroidal assumption the matrix B(9) becomes a block cyclic 
of the form, 

®11 ®12 ••• ® 1 M 

B(©) ^ 

®12 ®xi 

where each is a cyclic submatrix and depends on the 
parameters ^ For example, when the neighbour set of 


dependence 

N is 




, (e,x), (1,0), (0, 

-1)^ then we have 

hi "" 

: eye 



®12 = 

: eye 

( — 1 

h 

O 

• 

• 

• 

0 3 

®1M ~ 

! eye 


.. 0] 

and =s 

0 for 

jl ^ fOi. |M • 



To ensure stationarity, the coefficients ©j^ must obey 

i*s(a) 4 (1- ^ o^V-s feJi- (4.17), 

where, 

« col [exp[f-l ^ (s-l)^r, x ^n]-, I =(;,!) 
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In a similar manner, we can construct the simultaneous ARMA 
or SARMA model. 

Conditional Markov (CM ^ Model * 

In previous chapter we discussed the properties of 
a CM model defined over an infinite lattice. Again, we can 
have CM models defined over a finite lattice. Vie can 
also tackle the problem of edges by toroidal assumptions. 
The vector matrix relation in CM model is 

HC») I = e (4.13) 

where H(0) is a block cyclic and symmetric, AHere to 
ascertain stationarity we need 

•‘s’ 4 <**8^ > 0 V s 6-0- 

0 - col [0^, r^N] 

and 

cig = col [cos I^(s-.l)^r, r^Ng] 
where N is the asymmetrical half of N^i.e,, 

N = U fls 



Us = 0 
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4.7 PROPERTIES OF FL MODELS 

Considering the matrix B{©) defined in eqn. (4.X4) 
we can derive some informations about FL SAR model from its 
toroidal property. 


l) Since B(©) is a block cyclic matrix ^ hence the eigen- 
vectors P'jLj(O) of B(©) are the fourier vectors f.-» where 


^ (Mxl vector) 

and 

fj =s col [if^j# ] ( Mxl vector) 

with ss exp ^ ] 

So given the set of parameters 
the eqn. for as 


hi 


we can write 




= 1 + 


E 

(P,q)6N 




where k = (p>q) 


(4.19) 


2) Again from the vector-matrix relation we can write the 
covariance Q of Y is as 


Q cov [y] a [b'^(©) B(0) j:> 
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R(k,i) =: autocorrelation at lag {k,l) is given by 
R(k>0 


M . ft 

^>i Aj (4-20) 

The discrete spectral density as defined by the DFT C2D) of 
R(kj^) and then for such block cyclic matrix B(d) we have 


=f/ll 


(4.21) 


± 


llt^y(e)ll 


•2 >®'i = ¥ 

= 22e1 

^2 M 


CM Models: 


Since in FL toroidal CM model H(0) is a block-cyclic 
symmetric matrix, some proporties of CM(FL) models are 

1) Similar to FLSAR model, the matrix H(0) has fourier 
vectors as its eigenvectors* If we assume that ptj (©) is 
the eigenvalue of H(0) corresponding to f^j we have 

2) Without any loss of generality we can assume a = 0. 

Then, y » [H(»)r^ e^^ 


(4.22) 
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where e is a correlated sequence. 

To obtain the covariance of x we raultipily eqn, (4*22) by 
X* and take expectations on both sides. 


Q = [H(e)]-^ E[y'‘^e] 

-1 

= CH{e)] •>> 


(4.23) 

from eqn, (4.6). 


and covCe] » E [ee^] = H(9) 

(4.24) 


4*8 RELATIONSHIP BETWEEN SAR AND CM MODELS 

/ 

In general, two FLRF models defined over same 
lattice JX are said to be equivalent if their corresponding 
output vectors are identical. Unless otherwise stated^ 
equivalence means equivalence in all respects* 

Direct comparison of SAR and CM models are hampered 
by the fact that input e to the CM system is correlated. It 
would be desirable to modify CM model to a new model where 
the input is uncorrelated like that of SAR. This new 
representation is also helpful in synthesis of image 
following a given Cm model. 

The CM model as given by eqn. (4,22) has the following 

representation 

mo) [£l = ('‘■25) 
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ifHii) 


M 

Z 

if J= 




a 


3^3 


(4.26) 


where u is the normal zero mean IID sequence and ifHC©) is 
real and positive definite. 


Eqn* (4,25) again may be written by using egn^ 
(4,26) as 


M 


^ i,j«l ^ /^^iij*(«) 


(4.27) 


So e in eqn. (4,22) is related to u by the eqn. 


e ss Yh(©) u 


(4,28) 


Even though eqn, (4,25) has a form similar to that 
of a SAR model but, it is not necessarily a SAR model, A 
FL CM model represented as 


^^^1»^2^ [y(s)] - Y^e(s) (4.29) 

and it is equivalent to a SAR system represented by eqn. 

B(0*)x = V^H (4,30) 

if G(Zj^, 22) passes a factorization given in eqn, (4,3l) 

- S 1 |A(zj^fZ2)! 1^ 

> k A(Zjl,Z2)A(Zjl"^, ) (4.31) 

If G(2j^fZ2^ follow eqn, (4.3l) but there exist 
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polynomials 0 Ct|^(zj, 22 ) and ® constant k 

such that 

GCzijZg)! lai(2j^,Z2)l 1^ = k} |a2(z^,Z2)| (4,32) 

then the corresponding CM model has an equivalent SARMA 
model. If G(zj^^Z 2 ) factorizable then the CM 

model has no equivalent simultaneous model, 

4,9 PARAMETER ESTIMATION 

Suppose we have already specified the model type and 
the neighbour set. Let 0 - (0, P) be the vector of unknown 
parameters, where 0 is the vector of parameters and 
is the noise variance. The problem of estimating ^ from 
the given image vector y is primarily a statistical 
problem. This problem can be solved either by least 
square (LS) method or maximum likelihood (ML) method. In 
the case of LS estimate the estimates are given by 

0 = [E z(s) z'^Cs)]-^ (E z(s)y(s) (4.33) 

and 

P “ ^ ^ (y(s) - 0 z(s))^ (4.34) 

where 

z(s) = coa Cy{s+r), r6N] (4.35) 

It has been shown by Ord [37] that this method does not 
give consistent result for nonunilateral neighbourhood. 
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Lik elihood Est.ijna-te 

In case of Fl images, the mlE gives a consistent 
and efficient estimates of the parameters. For such MLE 
we have to find the log likelihood function of the output 
vector Y,. To obtain such an expression we impose the 
condition that input IID sequence is a zero-mean Gaussian . 
which gives a Gaussian output as the system is linear. 

If we represent the system by its tranfer function Aiz^,Z 2 ) 

or by its matrix from B(0) then likelihood function comes 
out to be 


p(2 1 a .f ) = j e=!p[- k I A(z^,Z2)y(s)]2] 

I s ^ -/V- 


where J is the Jacobian 
or. 


P(x|§l»'P) = In [det B(©)3 - |- In (2iip} 


- 1 (y(s)-0^z(s))^ 

J s ^ A. 

since det [B{9)] = II (1- ) 

wc have 


(4.36) 


(4.37) 


In P(l I a, p) = ^n(l- g)- ln(2iip) 


2 (y(s)- 0^ z(s))^ 

s^JV 


(4.38) 
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The ML estimates are obtained by maximising eqn* 
(4*35) 0 and -p , Since log likelihood function is 

a quadratic in 0, the estimation may be done by numerical 
optimization schemes like Newton-Raphson etc*^ which are 
computationally expensive 

Here we shall describe an iterative scheme for 
estimation which is an approximate version of the MLE, In 
such a scheme we approximate expressions like ln(l+a) by 
a - a /2« The corresponding log likelihood function 
denoted by J(0, becomes 

(4.39) 

(4.40) 

(4.41) 

(4.42) 

(4.43) 
algorithm 


J(E, p « 4- 0.5 0^R0 - In ( 211 ^) 


E (y(s) 

1 


©‘^z(s))2 


Where, 


Z Cg , m X 1 vector 




R 


a E ^ ^s matrix 




col [cos ^ ( (2.''l)^£) » ^ ^^3 


S = col [sin ((£*"ii )^£)> ^€^.^3 


's ' ■' M 


After defining the matrices and vectors the 
for approximate MLE is as follows: 
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?t+l =■ (R - ^ S)-l (y - U ) 
= ft - J - 

and 

ft = ^ (y(s)- 0 T2(s))2 

' W S£Jl t ^ 

where 

5 = s z(s) z(s)'f 

S4/1, 

(mxm matrix) 


U = E z(s)y(s) 
s (1^ JL.. 

(mxl vector) 


( 4 . 44 ) 


( 4 . 45 ) 


( 4 . 46 ) 


( 4 . 47 ) 


The initial value 0 is chosen as 0 = U, and here m 

is the dimension of §.. 


CM Model ; 

An estimate with good consistency and efficiency may 
be obtained in QA models by ML estimates. But it is not 
easy to derive an explicit expression fo.r the log likeli- 
hood function for such CM models because Jacobian of the 
transformation matrix is difficult to be evaluated. By 
assuming a toroidal structure and Gaussian correlated 
sequence e(s), the log-likelihood function may be given as 
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pCi| i = E I In (1-20^^ )- £ InCaTt-^) 

- ly Z^H(©)y 

(4.48) 

We see that the exponential term in In p(^ ©,')?') is linear 
in 0 unlike that of SAR model, A computationally efficient 
and consistent estimation procedure is as follows 

t == C s q(s)q(s)^]~^ ( z q(s)y(s)) (4.49) 

and 

y"" = ^ /r ^ q(s))^ (4,50) 

where JX| is the interior region. This procedure is 
suggested by Woods [36], * In our study, estimation 

of Parameters for both SAR and CM followed the above 
raontioned scheme and experimental results show that the 
estimates are quite close to the actual values and moreover 
this scheme works well for non Gaussian cases also (see 
section 4,11) * 

4.10 CHOICE OF NEIGHBOURHOOD FOR SAR AND CM MODELS 

From ID time series modelling, it is known that 

order 

only a model of appropiate £ gives the best result for a 
set of discrete data for forecasting and prediction. A 
similar situation is true in 2D spatial processes, but, the 
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problem with 2D data set is the variety of models available. 
So if we want to reconstruct any image from the parameters 
extracted from it, then the quality of reconstructed image 
depends on how the underlying model is similar to the true 
model and here comes the importance of the choice of 
appropiate neighbourhood. 

Though, one way of picking up the appropiate model 
is by visual inspection, we give below a quantitative 
measure for this. One way is to do pairwise hypothesis 
testing. But such acheme is not consistent and also not 
transitive. If a model is preferred to C2 and is 
preferred to then it does not follow is preferred to 
• Such a decision rule is not consistent in the 

sense that the probability of choosing an incorrect model 
does not go to zero if we take infinite number of 
observations. The other choice is based on Bayesian 
decision rule [35], 


Suppose that we have three sets 1^2 end of 
neighbours having *^3 respectively. 

For each model we can assign one SAR model Cj^ 

y(s) = £ ©j^y(str) +1^0^ 

3C N|^ ' 

y(s) ~ I \ yj(s+r)+irp j^u(s) s 
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where the significance of such equations are already 
defined. 

The models i = 1,2,3, are mutually exclusive. 

By Bayesian rule we can say that the optimal decision rule 
lor minimising the average probability of decision error 
is to choose the model that maximises the posterior 
probability , where y . Is the observation vector* 

The quantity P(Cj^ly) is computed by Bayesian rule, 

■“ PCl|C^i)P(Cj^)/P(y). We can set P(C^) same for 
ail i in the absence of any contrary information, so that 


p(x jCj^) =/pCyl©ij^) p{0,^ |Cj^)ld0l d|> 

This expression, with the Gaussian assumption is similar to 
the expression which is obtained in the section on para- 
meter estimation* By assuming that p(0,'p/C^) is a-priori 
density then by using asymptotic information it can be 
shown that 


P( - m2/2 in H- i E 0^3 S^) 


ln{M^)+ in 


In general, computation of such an expression for all values 
of k is not computationally efficient so here we will give 
an approximate test statistic by dropping the terms due to 
prior densities. Then the decision rule for the choice of 
appropriate neighbour set is following. 
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Choose the neighbour set N if 

k 


min Y C, 7 I 

k { 


( 4 . 51 ) 


where 


'k 5^)+ M^lnf 


2 X 

+ In(Nr) ^ 




Qt,c ~ s, ^ s, + c c. ^ 

ks ks ks ks ks 


Ss [cos ((s-lj'^r), r£N] 

\s = [sin |2E ((s-l)'^r), r^N] 

Bayes Rule for Choosing CM Model ; 


Suppose, WG have three sets Nj^#N2 and of neigh- 
bours having 2m2, and 2m^ members (symmetric) respe- 

ctively- 


For any model Nj^, we can write. 


y(s) E y(s+r) + e(s), s 

” X ®lfr ^ ic s 

t£u. "• ] 

where the symbols are defined previously, and e(s) is 
correlated Gaussian sequence. Then in the similar way, 
of SAFtf wo can say that the approximate test statistics 


will be; 
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Choos 0 the set N if 

Ic 


ram 

k 


where. 


A 




«'ks)+ 


+ nij^ In (M^) 


(4.52) 


where , 

^ = col [0 , 1 c£n ] 

\ 

and 


«'ks “ r r^N, ] 


where is the asymmetrical half of N, , i.e., 

''k k 


^ "“k ®k 


•r N )- , N n ^ 

s,. \ Sj^' Sj^ 


’k 


In the next section where we will discuss the simulation 
recultf. , we will show how test statistics detects the 
actual Jiiodel for synthetic images for both SAR and CM models, 

4.11 EXPERIMENTAL RESULTS 


To see how the SAR and CM models and associated 
estimation scheme work, different images, using different 
models have been generated and displayed for visual tests. 


'X, 
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The fast schemes for generating toroidal synthetic images 
by a known SAR and cal models are given below. 

Generation scheme for MxM iimaf|e (S^) ; 

i) Generate a sequence of IID numbers (Gaussian or 
uniform )of length and arrange them to form a M 5 cM 
array m Let that array be U . 

ii) Given a set of neighbours q. find the 

r 

eigenvalue matrix of the system as 

1 + E 0 2.^ where r=(kr^ 

Aj r 6N ^ ^ ^ 

iii) Take the 2D D.F.T* of the array U. Let the 
transformed array be 

iv) Generate an array where 

Fy(i.J)= (Fu (4.53) 

v) The final image obeying the given SAR model is 

y = f"’^(F ) 

« ^ F^(i,j)/ix(i,j))jp3 

In actual practice we did all the transforms by the MSL 
routine FFT2 and FFRDR2* 

The above algorithm for generating image was used 
with different neighbourhood (pre-specif ied) and with IID 
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sequences with different densities. Parameters were 
estimated from the generated image by using approximate 
ML estimates, in any random number generation, we require 
some initial values (called ISEED in IMSL) and it has been 
found by experiments that estimated parameters change wi-Ui 
different ISEED value and that change is also quite 
random. 

Generation Scheme for um Image ((M) t 

Since the generating eqn. in CM model can be 
written as 

fH(0) X = 

hence the generating strategy for CM is same as that of SAR 
except the cjqn. (4,53) where in CM, eqn, (4,53) is modified 
to 

Fy -(Fu (4.54) 

Hence also we take uniform as well as Gaussian density 
for the IID sequence . 

Results for SAR and GM model parameter estimations 
lor differi'nt ISEED value will be given in tabular form. 

In the? next set of tables we will give the estimated 
pa’ranu,*ti.?rs witii the test statistics Cj^ for finding the 
optimal neighbourhood. Here we have considered only 
noncnusal neighbourhood with 2,4*6 and 8 parameters for 
both SAR and CM models. 



ACTUAL 


GAUSSIAN SAR ESTIMATED 
(6 Point) 


100 



1 • 1 111 

0.9982 

1.035530 

1.081630 


0.28 

0.2988 

0.3012 

0.2877 

®10 

0-12 

0.1372 

0.1464 

0.1220 

®«-l— 1 

■~0 .14 

0.1612 

“0.1589 

-0,1729 

‘^11 

-0.14 

“0.1612 -0.1589 

-0.1729 

®10 

0.12 

0.1372 

0.1464 

0.1220 

^l*"! 

0.28 

0.2988 

0.3012 

0.2877 









GAUSSIAN SAR (4 POINT ) 



=1.1111 

1.02123 

1.07213 

1.00213 



^10 

0.36 

0.4123 
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UNIFORM (SAR) 6 POINT 


== 1,1111 

1.021346 

0.998312 

1.021672 

0.28 

0.3087 

0.3054 

0.3124 


CM 

• 

O 

0,1423 

0.1482 

0.1480 e_io 

•**0 #14 

-0.1642 

-0.1601 

-0.1650 

-0.14 

-0,1642 

-0.1601 

-0.1650 

0,12 

0.1423 

0.1482 

0.1430 ©i«i ®11 

0.2B 

0.3087 

0.3854 

0.3124 

CM (6 POINT) GAUSS 

1,1111 1.2312 

L. 01245 

0.99824 


®-ii 

.2800 

0.2912 

0.2992 

0.3012 

1 

d> 

®-10 

0.1200 

0.1421 

0.1291 

0.1282 

®0-l ®01 

®-l-l 

-0.1400 

-0.1232 

-0.1501 

-0.1542 

^10 


-0,1400 

—0 # 1232 

-0.1501 

-0.1542 
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®10 

0.1200 

0.1421 

0.1291 

0.1282 


O 

'^1-1 

0.2800 

0.2912 

0.2992 

0.3012 
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11 


CM (4 P0INT)GAUSS 

actual estimated 


0.28 

0.2942 

0.2852 

0.2915 
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-0.1493 

-0,1413 
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-0.1493 

-0.1413 

0.28 
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1.1111 

1 .0243 

1.0823 
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CM (4 POINT ) UNIFORM 

0.28 
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0.2901 

0.2713 
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-0.1482 

-0.1391 

-0.1542 

-0,14 

-0.1482 

-0.1391 

-0.1542 

0.28 

0.2752 

0.2901 

0.2713 ( 

1.1111 

1.1623 

1.1192 

1.2109 


© 


- 1-1 


© 
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- 1-1 


1-1 


© 


-11 


© 


11 


© 


-11 


0 


11 
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Computation of Parameter for an Image (following SAR 
Model) when the real neighbourhood is N = (0^1), 

(1,0), (-1,0). Original ®oi^O,-l~ 


l.Illl 


Positions 



1 .0052 


1 -0686 


1.0143 


0.4005 

-0.1152 

-0.1152 

0.4005 


0.5616x10“^ 

0.4041 

0.5550x10”^ 

0.6550x10“^ 

0 .4041 

0.56l6xl0~^ 


0.6963xl0~^ 

0.3966 

0.7091x10"'^ 

-0.1067 

-0.1067 

0.7091x10“^ 

0.3966 


80.76 


587,48 


U,0) 


(- 1 , 0 ) 


( 0 ,- 1 ) 


( 0 , 1 ) 


- ( 1 , 0 ) 


(- 1 ,- 1 ) (<* 1 , 0 ) (- 1 , 1 ) 



( 1 ,- 1 ) ( 1 , 0 ) ( 1 , 1 ) 


(- 1 ,- 1 ) (- 1 , 0 ) (- 1 , 1 ) 


230.07 (0,-1) 


( 0 , 1 ) 


( 1 ,- 1 ) ( 1 , 0 ) ( 1 , 1 ) 


0.6963x10 
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Confutation of the parameter Cj^ for an image (following CM 
model) when the real neighbourhood is N= (1,1), (-1,-1), 
(-1,1). Original “0*1^ = 1.1111 

®11 ®-l-l = 0*28 



0 a 


Positions 


1.1292 

©11=0.3221 

420.72 

-1-1 

11 

1.4213 

*-11=0.2121 

621.63 

1-1 

-11 

1.0321 

©_.11=0.1512 

©11 =0.2922 

121.43 

-1-1 

-11 



1-1 

11 


0^q=O.O212 -1^ 

1.1021 412.76 

9(^1=0.0092 Q^-L 01 


10 
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Real Image Synthesis using SAR Models ; 

Our discussions so far have been concerned with 
synthetic image patterns generated by known SAR models. 

We conclude this chapter by studying the appro piateness 
of the SAR models to real images. For this purpose, we have 
taken a 32x32 portion of Monalisa image. Since the entire 
image may not be homogeneous we have divided the entire 
image into 4 blocks of size 16x16. Parameters (©,^) 
were estimated for each block by the approximate ML estimate 
scheme with 12 point neighbourhood. The neighbourhood 
includes the points (1,0), (1,1), (0,1), (-1,1), (-1,0), 

(-1,-1), (0,-1), (1,-1), (2,0), (0,2), (-2,0), (0,-2). We 
have also estimated the parameters considering 8 point 
neighbourhoods. In that case, the neighbourhood N contains 
( 1 , 0 ), ( 1 , 1 ), ( 0 , 1 ), (- 1 , 1 ), (- 1 , 0 ), (- 1 ,- 1 ), ( 0 ,- 1 ), ( 1 ,- 1 ). 

In the case of 8 point neighbourhood the quality of the synthesiz! 
synthesized image was not good. From the estimated para- 
meters and sample means, we compute the residuals u as 


w 



[H(0) 1 - a 1] 


where x block image and a is the sample mean. From 

this sequence w we have generated a set of random numbers 


whose histogram was approximately same as the histograms 
of the residuals. The new sequence and the parameter 
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set were then used to synthesize the image* The synthesized 
image closely resembles the actual one except the edge 
effect of the block operation, 

4.12 DISCUSSION 

In this chapter we have given a complete methodology 
for analysing a given set of spatial data. We have also 
discussed the estimation scheme for SAR and CM models from 
the computation of test statistics. We have not said 
any thing about the superiority of SAR over CM model or 
vice-versa as the appropiateness of one class of model 
is entirely dependent on the data itself. 




CHAPTER 5 


TEXTURE SEGMENTATION BY SAR MODEL 

5.1 PRINCIPLE OF SEGMENTATION 

Segmentation is a process that breaks up a con^site 
scene into its constituents. Though there are several 
approaches for tackling this problem [3 ] y®'^ 1® an 

active field of research because of its importance as the 
first processing step in any practical scene analysis 
problem* In natural composite textures (e.g. combination 
of land and sea) we may or may not have any sharp boundary 
as the line of demarcation. In such cases where the 
contour is missing^ it has been observed that statistical 
measure gives an estimated contour which is very close to 
the original one. Here cementation by gray level 
thresholding will not work satisfactorily because there 
may be two regions having the same brightness average 
but slightly different variances so that it is not distin- 
guishable. In a composite image (texture) having two 
homogeneous regions we assume that all the pixels in one 
region are having same statistical properties. In fact, a 
composite texture may have more than two homogeneous 
regions but solving a two region problem can always be 
generalized for n-region textures. 
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Though there are several statistical segmentation 
schemes, the proposed scheme for segmenting a random field 
composite texture based on spatial interaction model is 
superior to them in the sense that it is computationally 
very simple and to some extent accurate. In this chapter 
we will discuss the principle of such segmentation scheme 
and the associated results, 

5,2 APPLICATION OF RF NIODEL FOR SEGMENTATION 

Given a composite texture X, we assume that for each 
homogeneous part of it we can have a RF model and model 
parameters are different from one region to another. Let 


X = U X^ (5,1) 

where X^^ and X 2 are the homogeneous random textures. Since 
and X^ are following a toroidal SAR model, we write 


y]:(3) 






= ■^^u(s).s(L'i-j 


(5,2) 


y^(s) + Z y5L(s+k^)=YfiU(s) sl-Aj 


k^^Ni *"1 


imilarly, 


X^s y2(s)+ S ©j^ y2(s+k2)=ff2'^(s) 




X. 
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where the various parameters appearing in eqn, (5.2) and 
eqn. (5.3) are explained in the previous chapter. The 
entire image is composed of two regions} region-1 and 

'f 

After estimating 
the parameters for each region we use them to find the 
similarity of them inside a region. The segmentation 
procedure begins by breaking the entire image into smaller 
Square blocks. We then consider each block image as an 
independent one and the parameters for each block is esti- 
mated, The same procedure is followed for all such blocks. 
We thus obtain a set of parameter for each block. The 
final step is to cluster these set of data by any standard 
scheme. The clustering of these data will determine the 
segmenting line. Though the entire image may be correlated, 
for the purpose of segmentation, we assume the block 
images to be independent. This is not a valid assumption 
and will give some error in the final result. But 
experience has shown that this error is negligible with 
moderate block size (8x8) and it varies inversely with the 
window size. The entire process can be visualized by 
Fig, 5,1* 

Hero we consider a MxM image with windows of size 
NxN* So, effectively, we will have K such blocks where 


respecti- 


region-2, represented by f 1 ^ 

vely, where and ©j^V 



XIO 


segmenting line 

Left portion of line 
with and 

right with (© 2 » (^ 2 ^ 


Fig. 5,1: A Block Composite Image 

K = M /N . Vi/e estimate the parameters of each block, 
treating them as independent, and arrive at K. sets of ^6, 

Then by clustering these parameters we can separate the 
non-uniform blocks. Since clustering in 2 region problem 
is a binary decision making process, the estimated curve 
will always be a staircase one. Evidently, the estimated 
curve will be close to the actual one if the block size 
is reduced but, we will show that it has other implications 
also. This method is tested with different segmenting curves 
like straight lines, parabolas , circles and for each case 
satisfactory results are obtained, 

5.3 PRINCIPLE OF CLUSTERING 

In this section we will discuss the principle of . 
clustering and the algorithm used. Clustering is also 
known as unsupervized learning. Roughly speaking, clustering 
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procedures yield a data description in terms of clusters 
or groups of data points that posses strong internal 
similarities. There are several methods of clustering for 
a set of data but all are heuristic in the sense that they 
do not possess any mathematical properties but are easy 
to use [38 3 < Since the clustering procedures have no 
mathematical foundations, they depend on the nature of the 
data. For all algorithms at least some properties of the 
data should be known a priori. 


In this study, the algorithm used for the segmentation 
procedure, is known as K~mean algorithm [39]. In this 
algorithm we assume that the number of clusters at final 
stage is K. This is one approach of general iterative 
optimization scheme. We will describe the algorithm and 
also give an example of it. 


The K-mean algorithm is based on the minimization 
of performance index which is defined as the sum of square 
errors. This algorithm consists of the following steps. 


1* Choose K initial cluster centres z^(l), Z 2 (l)t«« •jZj^.Cl) 
These are arbitrary and are usually selected as the first 
K samples of the given sample set, 

< 1 


2. At the Kth iterative step, distribute the samples 
airong the K cluster domain using the relation 



« 
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X^Sj(K) if 

llX - z.{K) n < j j X - z^(K)lj, Vi^l,2,...,K 

MK (5.4) 

where Sj(I<.) denotes the set of samples whose cluster centre 
is Zj(K).Ties in expression (5.4) are solved arbitrarily. 

3. Compute the new cluster centre Z^(K+l), V" 

•J 

from the result of step 2 such that the sum of square 
distance from all points in S.(K) to the new cluster centre 
is minimised. In other words^ the new cluster centre 
Zj(K+l) is computed so as to minimised the performance index 
Jj, where, 

J-l A S 1 |X-Z-(K+l)i 1^ j=l,2,...K 

^ ^xes-(K) ^ 

^ (5.5) 

The Z.(K+l) which minimises J. is simply the sample mean 

J 

Sj(K.). So, the new cluster centre is given by 

2^(K+1) = ^ S , X , j=l,2,...,K (5.6) 

^ X S^.(K) 

where N. is the number of samples in S.(K). The name 
*K«mean* is obviously derived from the manner in which 
the centres are sequentially updated. 

4. If Zj(K.+l) = 4 j(K) for j = 1,2>...,K, 

the algorithm has converged and the procedure has 
terminated. Otherwise go to step 2. 



113 


Convergence of such a scheme is not proved and 
convergence depends oh the value of K and the type of data, 

V^e will also give an example of such scheme for 
a 2D data • 


Example ; 

1. Let K=2 and let 
us take 

Zj^Ci) = [^]=Xj^ 
and 

Z^(l) = [ ^ ]=X, 
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Xj^-Z^(l)l I<l IXjL-Zpd)! i Fig. 5.2s 2D Clustering 
and 

llXj-z^d)! 1<1 Ixg-z^d)! I 

So, 

Sji^(l) = ^^1*^3^ 

S^ ( 1 ) ~ ^ * ^4 f * * * ^20 ^ 

Update the cluster centres 

0.00 
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Zo(2) - ^ ^ 18 (X2+X4+..*+X2 o)- ^5*63 ^ 
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4* Since Z.(2) 4^ Z.(l) V"4» so go to (2) 
J J ^ j 

2. With new cluster centres we get. 


l|Xg - Zj^{2)|l<||X^- Z2(2)1|,^:^= 1,2,. ...8 

I 

l!Xj^~ Z2(2)1 |<1 Ix^- Z^(2)|j,y£= 9,..., 20 


Si(2) 

82(2) 


Xi ... Xgj 

X 9 ... X^Q 


3. Update the cluster centre 


Zi(3) = t (X^ + ... + Xg) = [];;25] 


^2(3) = fe-CXg+ ... + 


4. Since Z- (3) 4 Z.(2), j=l,2 so go to step 2. Step 2 

'J J 

yields same result as in last iteration so Z2^(4)=Zj^(3) 

and ~ 29(3). Step 3 also yields the same results. 

Since Z.(4) = Z.(3), j=l,2 the algorithm has converged 

and it has yielded the cluster centres 

-t r 1 • 25 -1 A -7 r 1 
^1 ^ 1 . 13 ^ ^2 " ^ 7 . 33 ^' 


This algorithm is very sensitive to the initial 
cluster centres but generally the first two data in the 
data set are taken as the cluster centres for the first 
iteration. 
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Though the above example deals with clustering 
of 2 dimensional vector but in practice for image segmen- 
tation we need clustering of feature vectors of dimension 
n+1 (for n-point neighbourhood), 

5.4 SIMULATION RESULTS 


In this section we will discuss the simulation results 
obtained by the above scheme. Various composite textures 
were generated using different SAR models. V^e have taken 
only 64x64 image but various window sizes were used for 
segmentation purposes. An empirical relation between the 
window size and error is given where error is defined in 

the mean square sense. If 0^ represents the parameter 

A 

vector for a fixed window ca and 0^ be its estimate then 
the error is defined as 


4 


(a) 


OX, 


£ 






(5.7) 


where N is the length of vector 0, 

In the following, we present the method adopted for the 
sciginontat Ion scheme and the method of error estimation for 
the; synthetic composite texture. The generated image has 
two parts or regions as shown in the Fig. 5,3, 
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Fig, 5,3: Composite Texture 


This image was actually generated by superposition 
of two different images such that, in the composite image 
left portion of the line (at 45°) has parameter 
and the right portion is with parameter 
actual simulation 9^ = [0,12 0.12]^ and - 1*1111 and 
9^ = [0.36 0,36]^ and ^ 2 ~ ^*9999 was used. We partitioned 
the entire 64x64 image into blocks of sizes 8x8 and 4x4, In 
case of 8x8 blocks we get 64 such blocks and 256 in case 
of 4x4 blocks. In Fig, 5.3 the line passes through the 
blocks 8,15,22,29,36,43,50 and 57 and they are really the 
composite blocks. Parameters were estimated for each blocks 
and all the blocks on the left side of the line yields nearly 
same parameters as is to be expected and same is true for 
lh«‘ blocks at the right side. After parameter estimation. 



we have a set of 64 (for 8x8) vectors where each vector 
has three components when we estimate a 2 point problem. 
These vectors are called feature vectors and they are used 
for clustering. Since we have considered symmetric neigh- 
bourhood, the 3 point vectors reduces to 2 point because of 
the synmetry. These vectors then can be plotted on a 2D 
plane. Graphically, the separations are distinguishable. 

As mentioned earlier, all clustering schemes are heuristic 
and they have inherent limitations in separating data points 
In the present case also, we end up with a few misclassified 
samples but they are easily detectable. V/e show a segmen- 
tation result for the same texture for 2 point and 4 point 
parameter estimation in Figure 5.4, 
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Fig. 5.4: Segmenting line estimated by different Sp^ 
models. 
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In Fig. 5.5, we have given the actual and estimated 
segmenting curves. In order to study the error in esti- 
mated parameter with window size,- we have taken a pure 
homogeneous image having © = and =s laill. Then 

the parameters were estimated for the image by considering 
windows of size 4x4, 8x8, 16x16, 32x32 and 64x64. The mean 
square error ^ ^ was also computed. This procedure was 
done for different window size and different ISEED values. 
From all such different error values we have taken the 
maximum error for a particular window size. As mentioned 
earlier the errors also vary randomly with different ISEED 
values. The errors were then plotted against the windov; 
size. As expected, the error was minimum for window size 
64x64, In Table 5.1, we present the error (maximum) for 
corresponding window size. 

Table 5,1 

WINDOV^ ERROR (Max.) 

64 0.9098x10“^ 

32 0.2912x10”^ 

16 0.3556x10“^ 

8 0.9455x10“^ 

4 0.293829 

0.352612 


2 
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The error vs. window size plot as shown in Fig. 5,6^ 
that there may be an exponential relation between 
error and window size. This experiment was carried out for 
homogeneous images with different parameter sets, lie have 
repeated the experiment by using different image sizes like 
4bx48 etc. In all such cases, we have found that error and 
window sire are related by the equation 

€w (5,8) 

vdiere I 

C. •" mean square error 
« window size 

For the Imago with 0 = [0,12 0*12]^ and = 1.1111, if we 
estimate it as a 2 point problem then the eqn. (5,8) 
becomes 

1 AKiA ur**'i,297 

tw = “ 

But estimating parameters for a 4 point neighbourhood for 
Idio s.ime image gave 

t.,;= 1.46 vri •9''^ 

By changing either image size or the ISEED value it has 
beet! soon that A and B change. From simulation of different 
liyrithotic textures wc have seen that A and B depend on i) 
number of parameters, ii) image size and iii) the ISEED 
value. AS the ISEED value can be varied randomly hence 
A and B al *50 vary in a random fashion. For this reason 
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we h.wc 'liven different plots for this error* vs, window 
size in I'irj, 5*7 to show the variation of error with ISEED 
and vvin(bw i'-lze. Apparently, it is impossible to give any 
in .,1 theinatinnl expression for A and B in terms of the above 
inuntioncd variables. Vie have computed the error of such 
a r-c'.hunu' for two cases 



prom 1.1 u‘ 
obia i 
ihv error 
teimf- 0 1' 


Taking (iify) as the feature vector 
Inking g, as the feature vector 

plots it is scon that the additional information 

by inrluding 'Pin the feature vector does not change 
^ computed 

a|){)reciably* We have also ^ discrepency rn 

the window si2e, where discrepency is defined as 


fol low:;. • 


If straight 

Un«- r«‘prrsfnts the actual segmen 



Fig. 5,8 j Region of 
Confidence 


ling line?, then the hatched portion 
is di?i inr-d to bo the region of 

utjnl idnnci* aS shown in Fig. 5,8 

+ - rnn bo neneralized to other segmenting curves, 
;jamt‘ oonoopts can bo genoiaxo.^ 

i, U..- ..■v.iiiwl.ud =ogmonting line passes through this region 
that the soymontation procedure has converged- 
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perfectly* But in actual practice, the segmenting line 
may deviated out of this region* Vsle define discrepency 
aS the ratio of number of blocks outside the region to the 
total number of blocks. This discrepency follows (Fig. 5*9) 
an exponential relationship with the window size and it 
is reduced by considering larger neighbourhood. One more 
thing to be noticed here is that the discrepency-v^indow 
size relation is independent of the segmenting curve. All 
these experiments mentioned above were done for various 
types of segmenting curves like st, line, circles, st, line 
with double slope etc. and all the results for them 
were consistent. 

ie conclude our study of texture segmentation by 
discussing some of the disadvantages of the scheme 
developed here. 

For the purpose of elaboration, we consider the two- 
region problem again. Let us assume that the two parameters 



A 

and, ^21 ^ ^ 
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lower value and upper value 
oI t-lio OGlinia'te of 0^^ respectively. 

If WG select 0j^, 0^ and the window size such that 
the following conditions are satisfied^ 


i) 

1 ^IL ”* ^2L ! 

1 >> 0 

ii) 

1 ®1H “ ! 

1 » 0 

iii) 

1 “ ®2L 1 

j » 0 


where jX^Xl nieans the absolute difference between the vectors 
X and ^ , then the estimates of 0i|^ and ^ are well separated, 
and the clustering scheme will work perfectly. But suppose, 
for a particular 0^^ and ©g, the conditions i) and ii) are 
satisfied but instead of iii) we have 


^IH ^ 

then at the time of clustering, we may find intersections of 
such two clusters for some data points. Though we can 
vis sign these points to any of these clusters arbitrarily, 
if the* number of such points is large, then clustering 
scheme may give erroneous segmenting line. That is to 
say that if the differences between the corresponding 
elements in the vectors ©^ and ©2 is not greater than a 
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threshold value this scheme may fail to detect the actual 
cc'ymontinq line. This is because of the limitations of the 
estimation rule. If for some 9^ and ^ the scheme fails 
to detect the line properly we may say that statistically 
the texture is not a composite one and it may be treated as 
a homogeneous one. 




FIG. 5- 7 THE ERROR VERSUS WINDOW FOR 
DIFFERENT ISEEO VALUE. 



FI6.5.9 DISCREPENCY VERSUS WINDOW FQR 


A 64 X 64 SYNTHETIC IMAGE . 



CHAPTER 6 
CONCLUSION 


In our study of texture modelling, we have considered 
only the statistical modelling concept arvi the special type 
of statistical model, considered here, is known as ‘random 
field’ model. In the random field models, we have considered 
both causal and noncausal neighbourhoods. Procedures for 
implementation of such models have been given in details. 

We have cited examples of application of such models to 
texture segmentation and real image synthesis. A complete 
methodology has also been given for frequency domain 
characterization of textures. The application of 2D time 
series modelling to textural image has been studied for 
different Gaussian textures. Vie have also studied the 
applicability and limitations of pseudo random arrays in 
image modelling, 

A new method for image segmentation for a composite 
image by SAR models is investigated. The area of image 
modelling is a vast in scope and the present study only 
considers one specific class of models. We have not specific 
class of models. We have not even touched upon the entire 
field of syntactic image modelling. 
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Some future areas of investigation are cyclostationary 
models and time series models other than linear Gaussian, 
Though linear non-Gaussian time series models axe very 
useful for modelling of some real data, the time series 
models, having random parameters are very useful for real 
life image modelling. A large amount of research into 
these models is now being concentrated on the construction 
and application of computationally efficient algorithms 
to determine order and obtain estimates of the unknown 
parameters which have desirable statistical properties [40] • 
The development of 2D counterpart of the varying parameter 
time series model may show a new direction for image 
modelling. 
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